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Preface 


A few decades ago, the advent of high-speed electronic digital computers gave tremendous impetus to all 
numerical methods for solving engineering problems, and made it possible to solve with good accuracy 
many problems which previously could only be solved approximately. Finite difference methods, applied 
manually, gave way to finite element methods, which are still one of the most versatile and widely used, 
particularly in structural and solid mechanics. In thermofluids, methods of the finite volume type tend 


to be preferred. 


Slower to develop have been boundary element methods, based on boundary integral equations. Initial 
development was largely in the hands of mathematicians, as the underlying mathematics are relatively 
sophisticated. It was engineers, however, who turned boundary element methods into practically useful 


and powerful techniques. 


The purpose of this book is to serve as a deliberately simple introduction to boundary element methods 
applicable to a wide range of engineering problems. The mathematics are kept as simple as reasonably 
possible. Computer programs form an integral part of the boundary element approach and they are 
treated as such in the text. Several programs suitable for use on desktops or laptops are presented and 
described in detail and their uses are illustrated with the aid of a number of practical examples. Problems, 


with solutions, are provided at the ends of the chapters, for readers to solve for themselves. 


The programming language used in the main text is Fortran. Although it is somewhat unfashionable these 
days for general programming purposes, Fortran is still very widely used in engineering computation. 
Matlab versions of the programs are also provided in Appendices. Full listings of all the programs, both 


Fortran and Matlab, are available for download here. 


A prior knowledge of either Fortran or Matlab is desirable. The level of continuum mechanics, numerical 
analysis, matrix algebra, vector analysis and other mathematics employed is that normally taught in 
undergraduate engineering courses. The book is therefore suitable for engineering undergraduates and 
other students at an equivalent level. Postgraduates and practising engineers may also find it useful if 


they are comparatively new to boundary element methods. 

The book is presented in two Parts. Part I started with a brief review of the problems encountered in 
engineering, showing that they of two broad types. It then described boundary element treatments of 
problems of the potential type, using both constant and quadratic boundary elements. This Part II is 
concerned with elastic stress analysis problems of the plane strain and plane stress types. 


Imperial College London Professor Roger Fenner 
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Notation 


The mathematical symbols commonly used in the main text are defined in the following list. In some 
cases particular symbols have more than one meaning in different parts of the book, although this should 


not cause any serious ambiguity. 


A area of a solution domain 

[A] square matrix 

Aj coefficient of matrix [A] 

a4, A2, A3 constants in general boundary condition Equation 1.83 
[B] square matrix 

By coefficient of matrix [B] 

[b] column vector of known coefficients 

C torsional couple 

C, specific heat 

Crxs Cy Cyxs Cyy free term constants 

D flexural rigidity of a flat plate 

E Young's modulus 

e strain 

Fp, Fp drag and pressure flow shape factors for downstream flow 


fi function of position in Poisson’s equation 

ho function of position in biharmonic equation 

G shear modulus 

g acceleration due to gravity 

g heat generated per unit volume 

H height of a channel or solution domain in general 

H total rate of heat conduction 

h surface heat transfer coefficient 

i nodal point number 

unit vector in the x co-ordinate direction 

Jacobian of transformation (from global to intrinsic co-ordinates) 
J nodal point number 

unit vector in the y co-ordinate direction 
ratio of outer to inner radius for a cylinder 
thermal conductivity 

nodal point number 

unit vector in the z co-ordinate direction 
maximum dimension of the solution domain 


element number 


i 

k 

k 

k 

M total number of boundary elements in a mesh 
m 

N total number of nodes in a mesh 

N, 


shape function 


ay 


direction of outward normal to the boundary of a solution domain 


outward normal vector to boundary 
Download free eBooks at bookboon.com 


J 
J 
K 
L 
n 
n 


3 

. 
3 
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components of n in the x and y directions 


=) 


unit outward normal vector to boundary 


> 


components of 7 in the x and y directions 


> 
tad 
< 


source/force point on a boundary 


Cia) 


pressure gradient in the zco-ordinate direction 
pressure 

source/force point 

volumetric flow rate 

field point on a boundary 

field point 

a function of intrinsic co-ordinate € 

radial co-ordinate 

distance between a source/force point and a field point 
vector distance between points 

unitvector in the direction between points 

boundary of a domain 

distance along a boundary 

vector tangent to a boundary 

ratio between lengths of successive elements on a boundary segment 
part of a boundary forming element m 
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3 


distance along a solution domain boundary 

temperature 

traction kernel function 

remote temperature of surroundings in thermal convection 
time 

traction 


Se es ee 


displacement kernel function 
displacement or velocity in thexdirection 


ee 
Ba 


displacement in the radial direction 


< 
cy 


displacement in the@direction 

mean value of u 

velocity component of a boundary in the z co-ordinate direction 
displacement or velocity in theydirection 

mean value of v 

width of a channel or solution domain in general 


TE TPUusw sé 


displacement or velocity in thezdirection 
global Cartesian co-ordinates 


P< >< 
IS 
Ni 


components of body forces per unit volume in the Cartesian co-ordinate 
directions 

X,V,Z Cartesian co-ordinates 

[x] column vector of unknown quantities 


a coefficient of linear thermal expansion 
a coefficient in mixed boundary condition, Equation 2.32 
B coefficient in mixed boundary condition, Equation 2.32 
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Subscripts 
c 

e 

i,j,k 

L 

m 

n 

PI 


T 
XV, zZ 
0 

1,2 
1, 2,3 


Superscripts 
* 


* 


an angle 

change in temperature 

a small distance 

local intrinsic co-ordinate within an element 
angle of rotation per unit length of a bar in torsion 
angular co-ordinate 

an angle 

permeability of a porous medium 

viscosity 

Poisson's ratio 

local intrinsic co-ordinate within an element 
dimensionless pressure gradient 
dimensionless flow rate 

density 

stress 

velocity potential 

fundamental solution for potential 

stress function 

stream function 

potential 

grad operator 

harmonic (Laplacean) operator 

biharmonic operator 


counter for nodes within an element 

von Mises equivalent (stress) 

nodal point numbers 

solution to Laplace’s equation 

element number 

direction of the outward normal to a boundary 
particular integral solution satisfying Poisson’s equation 
radial direction in polar co-ordinates 

direction along a boundary 

segment 

total solution (Laplace plus particular integral) 
Cartesian co-ordinate directions 

angular direction in polar co-ordinates 

inner and outer of two concentric circular boundaries 
nodes of a quadratic element 


effective value under plane stress conditions 
modified quantity 
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Some Program Variable Names 


The Fortran computer program variable names widely used in the programs and main text are defined 


in alphabetical order in the following list. 


A 
All 


Coefficients of matrix [A] 


Matrix diagonal coefficient A,(potential problems) 


AITXX, AITXY, AITYX, AITYY 


AJ 
AK 


Coefficients at the diagonal of matrix (elastic problems) 
Matrix coefficient A, 


First kernel function contributing to the matrix [A] (potential problems) 


AKXX, AKXY, AKYX, AKYY 


ALPERP 
ALPERP2 
ALPHA 
ALPHAN 
ALPHASEG 
ALSEG 
ANG 
ANGFIR 
ANGSEG 
ANGSTORE 
AROW 
AROWX 
AROWY 
BDPSI 
BETA 
BETAN 
BETASEG 
BY 

BK 


First kernel functions contributing to the [A] matrix (elastic problems) 
Perpendicular distance from centre of curvature to mid point of a segment chord 
Square of ALPERP 

Element values of constants in mixed boundary conditions 

Nodal point values of constants in mixed boundary conditions (quadratic elements) 
Boundary segment values of constants in mixed boundary conditions 

Length of a segment chord measured between its end points 

Angular position of current end point on a curved boundary segment 

Angular position of first end point on a curved boundary segment 

Angle subtended at centre of curvature by a curved boundary segment 

Angular positions of end points on curved boundary segments 

Array storing element node contributions to the [A] matrix (potential problems) 
Array storing element node contributions to the [A] matrix (elastic problems) 
Array storing element node contributions to the [A] matrix (elastic problems) 
Coefficient of right hand side vector (matrix [B] times vector of knowns) 

Element values of constants in mixed boundary conditions 

Nodal point values of constants in mixed boundary conditions (quadratic elements) 
Boundary segment values of constants in mixed boundary conditions 

Matrix coefficient B, 


Second kernel function contributing to the [B] matrix (potential problems) 


BKXX, BKXY, BKYX, BKYY 


BK2 


Second kernel functions contributing to the [B] matrix (elastic problems) 
Non-singular part of second kernel function when P is the current element node 


(potential problems) 


BK2XX, BK2XY, BK2YX, BK2YY 


Non-singular parts of second kernel functions when P is the current element node 


(elastic problems) 
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BROW Array storing element node contributions to the [B] matrix (potential problems) 
BROWX Array storing element node contributions to the [B] matrix (elastic problems) 
BROWY Array storing element node contributions to the [B] matrix (elastic problems) 

BTX, BTY Coefficients of right hand side column vector (matrix [B] times the vector of knowns) 
CASE Alphanumeric plane stress or strain problem type 

D Perpendicular distance from P to the element containing node Q 

DPSI Nodal point values of the potential gradient solution to Laplace’s equation 

DPSIPI Nodal point values of the particular integral potential gradient function 

DPSIPIM Values of the particular integral potential gradient at the nodes of each element 
DPSISEG Values of potential gradient applied as boundary conditions to the boundary segments 
DPSISTORE Temporary store for potential gradient 

DPSIT Nodal point values of total potential gradient (Laplace plus particular integral) 
DRDN Rate of change of radius with distance along normal to boundary 

DUDZ Rate of change of displacement u with € along element 

DVDZ Rate of change of displacement v with € along element 

DZDE Jacobian of transformation from intrinsic co-ordinate € to y 

E Young’s modulus 

EGL Values of the intrinsic co-ordinate at the Gauss points (logarithmic quadrature) 
ELENGTH Lengths of the elements 

ESS Direct strain along boundary 
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ESTORE 
ETA 
EVX 
EVY 

Fl 
FLOWELEM 
FLOWIN 
FLOWOUT 
FLOWSEG 
FXELEM 
FXSEG 
FYELEM 
FYSEG 
HX 

HY 

I 

11, 12, 13 
IBC 
IBCD 
IBCE 
IBCM 
IBCN 
IBCP 
IBCPC 
IBCS 
IBCU 
IBOUND 
IC 
IDIRPC 
IEEND 
TEP1 
TEP2 
IFIRST 
IFLAG 
IGAUSS 
IINT 
ILAST 
IN 
IROW 


Stored value of Young’s modulus 

Intrinsic co-ordinate 4 

x component of the vector along an element 

y component of the vector along an element 

Constant function fin Poisson's equation 

Potential flow across an element 

Total potential flow into the domain 

Total potential flow out of the domain 

Flows of potential across the boundary segments 

Force on an element in x direction 

Total force on a boundary segment in x direction 

Force on an element in y direction 

Total force on a boundary segment in y direction 

Interval between points in the x direction used in domain integration 
Interval between points in the y direction used in domain integration 
Node counter 

Numbers of the three nodes of a quadratic element 

Type number of boundary conditions applied to the (constant) elements 
Counter for segments subject to applied potential gradient boundary conditions 
Type number of boundary conditions applied to the (quadratic) elements 
Counter for segments subject to applied mixed boundary conditions 
Type number of boundary conditions applied to the nodes (of quadratic elements) 
Counter for segments subject to applied potential boundary conditions 
Counter for point displacement constraints 

Counter for segments subject to applied stress boundary conditions 
Counter for segments subject to applied displacement boundary conditions 
Counter for boundaries 

Case number for logarithmic Gaussian quadrature 

Direction numbers of point displacement constraints 

Counter for element end points 

Counter for first end point of an element 

Counter for second end point of an element 

Numbers of first nodes on the segments 

Flag for ill-conditioning of the [A] matrix 

Counter for Gauss points 

Counter for internal points 

Numbers of last nodes on the segments 

Counter for nodes within an element 


Number of row in the [A] matrix 
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ISEG 
ISEGBC 
ISEGELEM 
ISEGEND 
ISEGMAX 
ISEGMIN 
ISEND 

IT 

IX 


MAXNB 
MAXNEL 
MAXNEQN 
MAXNNP 
MAXNPC 
MFIRST 
MLAST 
MMAX 
MMIN 
NBCD 
NBCM 
NBCP 
NBCPC 
NBCS 
NBCT 
NBCU 
NBOUND 
NEEND 
NEL 
NELB 


Segment counter 

Segment numbers for a particular type of boundary condition 
Segment numbers for elements 

Segment numbers for element end points 

Number of last segment on current boundary 

Number of first segment on current boundary 

Counter for boundary segment end points 

Number indicating type of Gaussian quadrature (normal or logarithmic) 
Counter for points in the x direction used in domain integration 
Maximum value of IX 

Counter for points in the y direction used in domain integration 
Maximum value of IY 

Node counter 

Jacobian of transformation from global to local intrinsic € co-ordinate 
Maximum number of columns in the extended [A] matrix 

Element counter 

Numbers of the elements adjacent to the first node of each element 
Numbers of the elements adjacent to the third node of each element 
Maximum dimension of the solution domain 

Maximum number of boundaries allowed by the array dimensions 
Maximum number of elements 

Maximum number of equations 

Maximum number of nodal points allowed by the array dimensions 
Maximum number of point displacement constraints 

Numbers of the first elements on the segments 

Numbers of the last elements on the segments 

Number of last element on current boundary 

Number of first element on current boundary 

Number of segments subject to applied potential gradient boundary conditions 
Number of segments subject to applied mixed boundary conditions 
Number of segments subject to applied potential boundary conditions 
Number of point displacement constraints 

Number of segments subject to applied stress boundary conditions 
Total number of segments subject to applied boundary conditions 
Number of segments subject to applied displacement boundary conditions 
Number of boundaries 

Number of element end points 


Number of elements 


Numbers of elements on the boundaries 
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NELSEG Number of elements on current boundary segment 

NEP1 Numbers of the first end points of the elements 

NEP2 Numbers of the second end points of the elements 

NEQN Number of equations 

NGAUSS Number of Gauss points 

NINT Number of internal points 

NNP Number of nodal points 

NNPB Numbers of nodal points on each of the boundaries 

NODE Numbers of the nodes of the elements 

NODEPC Numbers of nodes subjected to point displacement constraints 

NSEGB Numbers of boundary segments on each of the boundaries 

NSEGTOT Total number of boundary segments 

NU Poisson's ratio 

NX Number of internal points in the x direction used in domain integration 
NY Number of internal points in the y direction used in domain integration 
PI T 

PSI Nodal point values of the potential solution to Laplace’s equation 
PSIBOT Value of potential on the bottom edge of a rectangular domain 

PSHP Laplace equation potential at an internal point 

PSIIPT Total potential at an internal point (Laplace plus particular integral) 
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PSILEFT Value of potential on the left hand edge of a rectangular domain 
PSIPI Nodal point values of the particular integral potential function 


PSIRIGHT Value of potential on the right hand edge of a rectangular domain 


PSISEG Values of potential applied as boundary conditions to the boundary segments 
PSIT Nodal point values of total potential (Laplace plus particular integral) 

PSITOP Value of potential on the top edge of a rectangular domain 

PSIVAL Values of potential stored for domain integration 

Rl Distance from point P to the first end of element containing node Q 

RIX x component of radius vector from P to the first end of element containing Q 
RLY y component of radius vector from P to the first end of element containing Q 
R2 Distance from point P to the second end of element containing node Q 

R2X x component of radius vector from P to the second end of element containing Q 
R2Y y component of radius vector from P to the second end of element containing Q 
RATSEG Ratio between successive element lengths on current boundary segment 

RFN Value of function R(é) 

RSEG Radius of curvature of current boundary segment 

RX Component in x direction of unit radius vector from P to Q 

RY Component in y direction of unit radius vector from P to Q 

SD Shape function derivatives for quadratic elements (normal quadrature) 

SDL Shape function derivative values for quadratic elements (logarithmic quadrature) 
SF Shape function values for quadratic elements (normal quadrature) 

SFL Shape function values for quadratic elements (logarithmic quadrature) 

SEN Shape function value at a Gauss point 

SIGE Nodal point values of von Mises equivalent stress 

SIGNN Nodal point values of direct stress normal to boundary 


SIGNNSEG Boundary segment values of direct stress normal to boundary 


SIGSN Nodal point values of shear stress along boundary 
SIGSNSEG Boundary segment values of shear stress along boundary 
SIGSS Nodal point values of direct stress along boundary 

STORE Stored values in the boundary condition application process 
TITLE Alphanumeric title for the problem (maximum 80 characters) 
TMX Element nodal point values of traction in x direction 

TX Nodal point values of traction in x direction 

TMY Element nodal point values of traction in y direction 

TY Nodal point values of traction in y direction 

U Nodal point values of displacement in x direction 

UELEM Element values of displacement in x direction 

UNGX x component of the unit normal at a Gauss point 

UNGY y component of the unit normal at a Gauss point 

UNMX x components of the unit normals at the nodes of each element 
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UNMY y components of the unit normals at the nodes of each element 


UNX x components of the unit normals at the nodes 

UNY y components of the unit normals at the nodes 

USEG Boundary segment values of displacement in x direction 

UV Nodal point values of computed displacements (or tractions) 

Vv Nodal point values of displacement in y direction 

VELEM Element values of displacement in y direction 

VSEG Boundary segment values of displacement in y direction 

WG Values of the Gaussian weighting factors (normal quadrature) 

WGL Values of the Gaussian weighting factors (logarithmic quadrature) 
XC x co-ordinate of the origin for the particular integral function 
XCENT x co-ordinate of the centre of curvature of a curved boundary segment 
XEEND x co-ordinates of the element end points 

XFIRST x co-ordinate of first end point of current boundary segment 

XINT x co-ordinate of an internal point 

XLAST x co-ordinate of last end point of current boundary segment 

XMID x co-ordinate of the mid point between the ends of a curved segment 
XNODE x co-ordinates of the nodes 

XP x co-ordinate of point 

XPOINT Global x co-ordinate of an internal point 

XQ x co-ordinate of Gauss point 

XSEND x co-ordinates of the boundary segment end points 

XX x co-ordinate relative to the origin for the particular integral 

YC y co-ordinate of the origin for the particular integral function 
YCENT y co-ordinate of the centre of curvature of a curved boundary segment 
YEEND y co-ordinates of the element end points 

YFIRST y co-ordinate of first end point of current boundary segment 

YINT y co-ordinate of an internal point 

YINTGL Values of y direction integrals stored for domain integration 

YLAST y co-ordinate of last end point of current boundary segment 

YMID y co-ordinate of the mid point between the ends of a curved segment 
YNODE y co-ordinates of the nodes 

XP y co-ordinate of point P 

YPOINT Global co-ordinate of an internal point 

YQ y co-ordinate of Gauss point Q 

YSEND y co-ordinates of the boundary segment end points 

YY y co-ordinate relative to the origin for the particular integral 

ZETA Intrinsic co-ordinate & 

ZG Values of the intrinsic co-ordinate § at the Gauss points (normal quadrature) 
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5 Boundary Element Analysis of 
Plane Elastic Problems 


In this chapter a form of boundary element analysis for two-dimensional elastic stress analysis problems 
such as those outlined in Chapter 1 is presented. Only quadratic elements are considered. Three- 


dimensional problems are discussed briefly in Section 7.3. 


In Sections 1.2.5 and 1.2.6 it was shown that both plane strain and plane stress problems can be described 
in terms of a fourth-order biharmonic differential equation for Airy’s stress function. The relevant 
Equations are 1.71 and 1.77, both of which are special cases of Equation 1.85. Clearly, such problems are 


fundamentally different from potential problems, which are governed by the second-order Equation 1.84. 


Using an Airy stress function approach to solving plane elastic problems can be very appropriate 
when the boundary conditions are defined in terms of stresses. In more general problems of practical 
engineering interest, however, boundary conditions are typically defined in terms of a mixture of stresses 


and displacements, and a different approach is to be preferred. 
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Following the methodology developed in Chapter 2 for potential problems, what are required are a 
fundamental solution, equivalent to Equation 2.13, and a relationship equivalent to Green’s symmetric 


identity, expressed in Equation 2.25. 


5.1 Fundamental Solution 


In the case of potential problems, the fundamental solution (Section 2.1) was in practical terms that due 
to a source of potential concentrated at a point in a solution domain of infinite extent in all directions. 
A mathematical requirement of a fundamental solution is that its value is singular (goes to infinity) at 


a point — the point where the source is located. 


In elasticity problems the corresponding fundamental solution is that due to a force concentrated at 
a point in an infinite domain, which again has the property of singularity at the point concerned. 
Immediately there is a substantial difference between the two classes of problems. The fundamental 
solution for potential problems involves only a scalar quantity: the source strength, and the effects of the 
point source in terms of potential distribution will be the same in all directions moving away from the 
point. In contrast, the fundamental solution for elasticity problems involves a vector quantity: the point 
force has both magnitude and direction. The effect of the force in terms of stresses and displacements 


is certainly not the same in all directions moving away from the point. 


In three dimensions, the fundamental solution is truly that due to a concentrated point force. In two 
dimensions, however, the problem is either one of plane stress, when the solution domain is very thin 
in the third dimension normal to the domain (Section 1.2.6), or plane strain, when the solution domain 
is very thick in the third dimension (Section 1.2.5). In either case, the fundamental solution needs to be 
thought of as that due to a force uniformly distributed along a line through the domain thickness in the 
third dimension. The force is a line force per unit thickness, which is applicable to either plane stress or 


plane strain. In the plane of the solution domain it still appears as a force acting at a point. 


Rather than deriving the fundamental solution from first principles, the approach adopted here is to state 
the result, and then show that this satisfies all the relevant requirements. Figure 5.1 shows a plane with 
both Cartesian (x,y) co-ordinates and polar (r,@) co-ordinates, both with the same origin, and the 
angular co-ordinate measured anti-clockwise from the r direction. Stress components o,, (radial direct 
stress), 0,, (hoop direct stress) and o,9 (shear stress) are shown acting on a small region of the domain. 
As in Equation 1.1, shear stresses are complementary, so that 0,g = 09,, and 0;g is shown acting on 


all four faces of the region. 


Download free eBooks at bookboon.com 


Figure 5.1 Stress components in polar co-ordinates 


Due to a line force of unit magnitude (per unit length) acting in the x direction at the origin, the three 


stress components at distance r from the origin and angular location 6 are given by 


(34+v*) cos 6 


oy = — SEs (5.1) 
(1-v*) cos 0 

ogg =e (5.2) 

Op = ee (5.3) 


The displacement components, y,. in the radial direction and ug in the hoop direction, are given by 


_ (1+v*)(3-v*) cos 6 d 
Uy, = PGI In (5) (5.4) 
ug = =2latv'r-a+v9@B-vyIn(5)| (5.5) 


where the length d is arbitrary, and remains to be chosen. It appears because the problem (of a 


concentrated force in an infinite medium) has had no displacement boundary conditions defined. 


The parameters E* and y* in these equations are the effective Young’s modulus and Poisson’s ratio for 
the elastic material of the domain. Their definitions depend on whether the problem is one of plane 


stress or of plane strain. For plane stress 
E*=E and vr=v (5.6) 


while for plane strain 


_ E 
~ (1-v?) 


E* and (5.7) 
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In other words, under plane stress the effective properties are the actual properties, while under plane 
strain they are modified according to Equations 5.7. A plane strain problem can be treated as plane 


stress, provided the modified properties are used. 


Figure 5.2 Stresses acting on a circle about the point of force application 
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Figure 5.2 shows a circular region of radius r of the domain centred at the origin. The stresses acting 
on a small arc AB of the outer surface of the region of angular extent dO are o,,. and 0,9 as shown. 
For unit thickness of domain in the direction normal to the plane shown, the forces on AB in the radial 
and tangential directions are o,,. x rd@ and 0,9 X rd@, respectively, and the total force on the region 


in the negative x direction, which should be equal to the applied line force, is 


Force = or cos 0 +d, sin@)rd@ (5.8) 
- . [(3 + v*) cos? 6 + (1 — v*) sin? 6]d0 (5.9) 

Now 
i cos*@d@ = [a sin?6d0 = 1 (5.10) 


so that Equation 5.9 becomes 
1 ; ‘ 1 
Force = ,,Gtv)at+d-v*)t]=7GB+1)=1 (5.11) 


The force is of unit magnitude, as required. An integral expression similar to Equation 5.8 can be used to 
show that the total force in the y direction is zero. Note that the requirement that the total forces on the 


circular region be constant, independent of radius, means that all the stresses are proportional to r7!. 


In plane polar co-ordinates, the strains defined in terms of displacements (equivalent to Equations 1.2 
and 1.3) are 


Ou, (1+v*)(3—v*) cos 6 
e- = — = - 5.12 
TF or 4nE*r ( ) 
uy , 1dug (1+v*)? cos @ 
Cag un bu 5.13 
60 r r Or 4nE*r ( ) 
1 du, 0 (“2) 2(14+v*)(1—-v*) sin @ 
erg = - r—|—) = ———_—_ 5.14 
ré r 00 - or\r 4nE*r ( ) 


Identical results should be obtained by applying the elastic constitutive equations (equivalent to 
Equations 1.20 to 1.25) to the stresses defined in Equations 5.1 to 5.3. 


Plane stress 


Under plane stress conditions and in the absence of temperature changes, the direct stress o,, normal 


to the plane of the domain is zero, and 


os 0 


ee _c¢ (14+v)(3-v) cos 0 
err = 5 Orr — VI) =F ry 


4nEr 


[-3+v)-v(Q-v)] = (5.15) 


TEY 
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which, in view of Equations 5.6, is identical to Equation 5.12. Similarly 


cos 0 
AnEr 


= _ (1+v)? cos 6 
[(A-v) +v@+y)] = ee" (5.16) 


1 
€99 = 7 ee VO) = 


which is identical to Equation 5.13. Finally 


_ org _ 2(14+v) _ 2(1+v)(1-v) sin @ 
erg = =, So 


G E19 4nEr (5.17) 


which is identical to Equation 5.14. 


Plane strain 


Under plane strain conditions, and again in the absence of temperature changes, the direct strain e,, 


normal to the plane of the domain is zero, from which 


1 


ez, = E eee = V(Orr + 06 )] = 0 (5.18) 


O72 = V(G;r + 90 ) (5.19) 
Then, using Equations 5.7 


[O,- = Vv (096 + V0,, + V099 )] 


ty] 


16; - v(099 + zz) | — 


fy] 


Crr 


1 
= [(1 — v*)o,, —v(1 + v) a6] 


1 _ +v*)G-v") cos 6 


v 1 * 
= =|¢,, - —- 396 | = Flor — V"096] = ee (5.20) 
as in Equation 5.12. Also 
1 1 
C00 =F [F699 — VOzz + 0, )] = E [F969 — V(V0,, + VO9g + G;r)] 
1 2 
= E Ket = ogg — V1 + V) Or, | 
1 v 1 m (1+v*)* cos @ 
= 5 [400 — 35 or] = ge lee — Vo] = ey) 
as in Equation 5.13. Finally 
Or 2(1+7) 2(1+v*)(1—-v*) sin 6 
erg = s =~ Ore = a (5.22) 


as in Equation 5.14. 
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Equations 5.15 to 5.17 and 5.20 to 5.22 show that the expressions for strains, under either plane stress or 
plane strain conditions, obtained from the stresses via the constitutive equations, are identical to those 
obtained from the displacements. With the stresses also satisfying equilibrium with the applied line force, 
this means that Equations 5.1 to 5.5 represent the true solution for a line force at a point in an infinite 


plane, in other words the fundamental solution for plane elastic problems. 


5.1.1 Displacement kernel functions 


y ” 


Vv 
. 
q 


p 


Figure 5.3 Displacements at the field point 


These results now need to be generalised for line forces applied in both the x and y directions, to give 
the displacements and tractions in the same pair of directions. The term traction will be explained 
shortly. Figure 5.3 shows the displacements at a field point q produced by a line force of unit magnitude 
in the x direction at force point p. Displacements y, and ug are given by Equations 5.4 and 5.5. The 


corresponding displacement in the x direction at point q is 


u =U; COSO — Ug sind 


7 Ca (-) 7 sneha +v*)?-(1+v*)3-v*)In (5) 
oe 


Ar E* T Ark * 


= any les In (5) —1+ cos? a| (5.23) 


It is convenient to choose the arbitrary length d such that 


_ (tv*)? [(3-v*), (1 2 
— 4n E* lees In (7) aes a| (5.24) 
Now cos@ = f, sind =f, (5.25) 
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where 7, and f, are the components in the x and y directions of the radius vector of unit length in the 


direction from point p to point q. Equation 5.24 can be written as 


(14+v*)* [@-v*) 1 aD 
U = Uy (p,q) = > [EO In (2) + FA |r(p. 4) # 0 (5.26) 


Ux (p,q) is the displacement kernel function defining the displacement in the x direction at q due to 


a unit line force in the x direction at p. 


Again from Figure 5.3, the displacement in the y direction at point p is 
v =u, sind + Ug sind 


_ (1 + v*)(3 — v*) cos @ sind " (“) cos 8 sin@ la ty)? - (14098 —-v9In (| 


4m E* 4m E* 
1 #\2 : 1 *\2 es 
v = U,y (p,q) = cr cosé sin 0 =O BA, (5.27) 


Ux (p,q) is the displacement kernel function defining the displacement in the y direction at q due to 


a unit line force in the x direction at p. 
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For a line force of unit magnitude in the y direction at point p (Figure 5.3), similar results can be 
obtained. The angular position of field point q relative to the line of action of the force is now not @ but 
- (= - a); so sin @ becomes — cos 9, and cos @ becomes sin @. Using Equation 5.24, the displacement 


in the direction of the applied force at q is 


v= Uy (og) = HF 9(°) +29] =H [EC eR] 2 


rp, q) #0 
Using Equation 5.27, the displacement in the direction at 2/2 anti-clockwise from the direction of the 
applied force at q is 


=, EV? 2 Bes 
u=——7 sin 8 cos 


(1+v*)? 
4n E* 


and u = Uy, (p,q) = PF (5.29) 
The four equations, Equations 5.26 to 5.29, define the four displacement kernels. They can be expressed 
elegantly in a single equation if tensor notation is introduced. For present purposes, however, it is more 
convenient to have them in their more explicit, if more lengthy, forms. The first subscript of U indicates 
the direction of the unit line force at p, while the second subscript indicates the direction of the resulting 


displacement at q. 


5.1.2 Traction kernel functions 


The displacement kernel functions are expressed in terms of displacements at the field point in the global 
co-ordinate directions. Similar results are required for stresses. For this purpose, the concept of tractions is 


introduced. A traction is the resultant stress, or force per unit area, on a surface in a particular direction. 


¥ 


Pt 


Figure 5.4 Tractions and stresses at a surface within a domain 
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Figure 5.4 shows both the Cartesian stress components and the equivalent tractions,t, and t, in the x and 
y co-ordinate directions, acting on a surface of a small triangular piece of material within a domain. The 
unit outward normal vector to the surface, 7, is inclined at angle 6 to the x direction. For equilibrium 


of forces in the x direction (bearing in mind that the domain is of unit thickness) 
ty, X (AC) = o,, X (AB) + oy, X (BC) 
ty = Ox, COSO + Oxy SIND = Oxy hy + Oxy Ny (5.30) 


where fi, and Ay are the components in the co-ordinate directions of the unit normal vector, and AC, 
AB and BC are the lengths of the sides of the small triangle. Similarly, for equilibrium of forces in the 


y direction 
ty X (AC) = ayy x (BC) + oy, X (AB) 


ty = dyy sind + O,, COSd =Oyy fy + Oxy fy (5.31) 


y yyy 


Figure 5.5 Stress components in polar co-ordinates 


Figure 5.5 shows the polar co-ordinate stress components at field point q due to a line force in the x 
direction at p. Also the equivalent tractions in the global co-ordinate directions,t, and t,, acting on a 
surface of a small triangular piece of domain whose outward normal is inclined at angle y to radius r, and 
angle § to the x direction. For equilibrium of forces on the right-angled triangular piece in the x direction 


t, X (AC) =0,; X (AB) cos 8 — 0,9 X (AB) sin@ — ogg X (BC) sin@ + 0,9 X (BC) cos 8 


ty = 0;, COSY COS O — O,g cosy sin# — dgg siny sin@ + 0,9 siny cos 8 
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With the stresses defined by Equations 5.1 to 5.3, this becomes 


_ G+tv*)cos*@cosy (1—v*)sin* 6 cosy 


7 Arr Anr 


1 
aa _ 4% * 2 
= [((1 — v*) + 2(1 + v*) cos? 6] cosy 
t, = = [((1 — v*) + 2(1 + v*) cos? o\— (5.32) 


where 7 = cos y- The traction kernel function for defining the traction in the x direction at q due to 
nm 


a unit line force in the x direction at p is 


ir 


1 * AVA A d 
T.(0,q) = —Zo [A —v*) +20 + v")R AI SE r(p,q) # 0 (5.33) 


Similarly, for equilibrium of forces in the y direction 
ty X (AC) =a,, x (AB) sin 6 + 0,9 x (AB) cos 6 + ogg X (BC) cos 6 + 0,9 x (BC) sin é 
ty = 0,, Cosy sin @ + 0,9 cosy cos @ + dgg siny cos @ + 0,g siny sin @ 


v*) 


=a [(3 + v*) — (1 —v*)] cos @ sin@ cosy + siny (5.34) 


d— 
4nr Anr 
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Now siny = sin(é6 — @) = sind cos 8 — cos 6 sind 


A 


= fiyh, — ih, (5.35) 


where fi, and fy are the components in the x and y directions of the unit outward normal to surface 
AC at point q. The kernel function for defining the traction in the y direction at q due to a unit line 


force in the x direction at p is 


1 
4nr 


*\A A d * AA a Zn 
Try (2,9) = [2a ave het = = 1) 075 = fyfi,)| (5.36) 


r(p,q) #0 


For a line force of unit magnitude in the y direction at point p (Figure 5.5), similar results can be 
obtained. The angular position of field point q relative to the line of action of the force is now not @ but 


ui , , 4 P me 
= ( = 0), so sin @ becomes — cos 6, and cos @ becomes sin @. Also,f, becomes fy and A, becomes 


y 
—f,. Using Equation 5.32, the traction in the direction of the applied force at q is 


ee oe _ 4y* * Cae dr 
ty = rege v*)+2(1+ v*) sin Aes (5.37) 


| Qa 


r 


and Tyy @,q) = -—[(1 —v*)+2(1+ v*)t,F, | r(p,q) #0 (5.38) 


Qa 


n 


Using Equation 5.34, the traction in the direction at 1/2 anti-clockwise from the direction of the 


applied force at q is 


1 : dr pe Le x 
—t, = -——|2a + v*) (-sin@) cose (1 —¥")(—f, sin? =, (— cos @))| 


1 gyi dr fos se a 
i, = -—|201 + v*)sin@ cos 0 — — (1 —v*)(A, sin 6 — fy, cos )| (5.39) 
1 e\a a dr VDD An 
and Ty (p,q) = ae UT, = (=) Gig = f.Ay)| (5.40) 


r(p,q) #0 


Again Equations 5.33, 5.36, 5.38 and 5.40 for the traction kernel functions can be generalised in a single 


formula using tensor notation. 


The displacement and traction kernel functions defined here for elastic stress analysis are equivalent to 
those defined in Equations 2.13 and 2.34 for potential problems. There is a considerable difference in 
complexity, reflecting the vector nature of a concentrated force against the scalar nature of a point source. 
But, there are important similarities in that both pairs are functions of either In Gt) (displacements 
and potential) or r~! (tractions and potential gradient). The techniques of dealing with such functions 


developed for potential problems will be applicable to elastic problems. 


Download free eBooks at bookboon.com 


5.2 Boundary Integral Equations 


The relationship equivalent to Green’s symmetric identity for potential problems required to develop the 
boundary integral equation for elastic stress analysis problems is Betti’s reciprocal theorem. Suppose that 
an elastic body is subject to a number of forces, F;,, at particular points and in particular directions on 
its surface, and that the corresponding displacements at the same points and in the same directions as 
the forces are u;,. Suppose also that the same body is separately subject to another independent set of 
forces, F;, at the same points and in the same directions, producing corresponding displacements u;. 


Betti’s theorem states that 
Lik Fe Ue = Ln Fe Ux (5.41) 


It is a form of virtual work principle: the total virtual work done by the first set of forces moving through 
the second set of displacements is equal to the virtual work done by the second set of forces moving 


through the first set of displacements. 


For present purposes, the first set of forces and displacements can be those associated with the particular 
problem to be solved, and the second set with the fundamental solution. Also, a force F, which must in 
reality act over a finite area, is equivalent to a traction, t, and the summation in Equation 5.41 becomes 


an integral over the boundary surface, S$ 
Si teuj dS = J teupds (5.42) 
Suppose that some point within the solution domain serves as the origin and point of application of 


the concentrated forces (in the two co-ordinate directions) for the fundamental solution, as shown in 


Figure 5.6 


small circle 
radius € 


Figure 5.6 A solution domain, including a small region of radius surrounding the point p 


Now, the fundamental solution is not valid at the point p itself. Consider therefore a small internal 
boundary, S,, of radius ¢ centred at p. Equation 5.42 can be written for the region between S and S,, 


which excludes p where both displacements and tractions are singular 


Sy tyme dS + J, thupdS = fo thuz dS + J ty, ds (5.43) 
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The next step is to consider what happens as the radius € shrinks to zero. On the circle, a displacement 


kernel u; (Equations 5.26 to 5.29) takes the general form of 
* 1 
ui = f:(@)In(=) + A) (5.44) 


where f; (@) is zero for U,y and Uy,. If € is small the value of traction t;,, is effectively constant over the 


boundary S,, and equal to its value at point p. Therefore 
21 1 21 
[ wuas= af p@m(Zjed+a@ | A@ede 
Ss 0 0 


1 
= t(p)e In (=) C1 + ee Ce 


where C, and (C, are constants. The second term obviously goes to zero as ¢ > 0. Similarly, 


lim,_s9 E In (-)| = 0, and consequently 
Jz, teuzdS = 0 (5.45) 
Also, a traction kernel t; must be such that 


Sc te dS =1 (5.46) 
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if the traction is in the direction of the unit concentrated force at p (and zero if it is at right angles to 


the force). Therefore, if the force is in the x direction 

J;, teu ds = u(p) (5.47) 
and as the boundary S, shrinks to zero size at the point p, Equation 5.43 becomes 

u(p) = Jy tpupdS — J. thug dS (5.48) 


The concentrated force in the x direction at p creates displacements and tractions in both the x and y 


directions at all points on the boundary of the domain, as defined by the kernel functions. Therefore 
wr) =| [ee Que. Q) + ty (QUey .Q]AS@ 
S 
— f, [UCQ)T xx (DV, Q) + VQ)T ry (P, Q)]AS(Q) (5.49) 
Similarly, if a unit line force is applied in the y direction at p 
v(p) = [_ [eCQDU ye, 0) + ty QW yy, Q)]AS@) 
S 
— f, [UCQTyx @, Q@) + V(Q)T,y (v, Q)]d5(Q) (5.50) 
The point isa field point on the boundary which moves along the boundary as the integrations proceed. 
The displacement and traction kernel functions are known, so provided the values of both the 
displacements and tractions are known at every point along the boundary, Equations 5.49 and 5.50 
provide a means of calculating the solution displacements at any point within the solution domain. 
They can be differentiated to find strains at p, and hence stresses via the elastic constitutive equations. 
In engineering practice, stresses and strains at points within an elastic solution domain are rarely required. 
In all but a few classes of problems the largest individual stress components and the most intense states 


of stress are located on a boundary. So for present purposes, data at internal points are not computed. 


Let point p be taken to the boundary at a point P, and again consider a small region of exclusion with 


boundary S, of radius centred at P. While Equation 5.45 still holds, Equation 5.46 becomes 


Js. t; dS = constant (5.51) 
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if the traction is in the direction of the concentrated force at P, where the constant is no longer unity. 
Indeed, if the boundary at P is smooth, the value of the constant is %, because exactly half of the 
complete circle is within the solution domain. If the boundary is not smooth, however, the value of the 
constant is difficult to evaluate directly. Similarly, if the traction is in the direction at right angles to the 
direction of the concentrated force at P, Equation 5.51 also applies, where the constant is no longer 


zero. Equations 5.49 and 5.50 become 
Cox PUP) + Cry PIP) = | [Ea(QWex Ps @) + ty Dey P, Q)]AS@) 
— f, [UCQ)T ex (P,Q) + V(Q)T xy (P, Q)]d5(Q) (5.52) 
Cx (PIUCP) + yy (PYOCP) = J [te (QU yx P,Q) + ty(QUyy (P, Q)]AS(Q) 
— f, [uCQ)Tyx (P,Q) + v(Q)Tyy (P, Q)]d5(Q) (5.53) 


where Cy (P),Cxy (P), Cy, (P) and Cy (P) are constants. The free terms containing these constants 
contribute only to the coefficients at the diagonal of the relevant matrix in the numerical implementation 


of the method, and these coefficients can always be evaluated indirectly. 


Equations 5.52 and 5.53 are now a pair of boundary integral equations for point P as the point at which 
line forces of the fundamental solution are applied, the first for force in the x direction and the second 
for force in the y direction. They can be used to determine the distributions of the displacements and 


tractions over the boundary. 


5.3 Discretisation of the Boundary Integral Equations 


In general, Equations 5.52 and 5.53 cannot be solved analytically, and some form of numerical method 
must be employed. Boundary displacements and tractions are found not as a continuous algebraic 
functions of position along the boundary, but as numerical values at a finite number of discrete points 
on the boundary. The boundary may be subdivided into small pieces or boundary elements. Associated 
with each element are one or more of these points: the nodes or nodal points. The distributions of 
displacements and tractions over the elements are defined in terms of nodal point values by suitable 


interpolation functions. 
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Boundary integral Equations 5.52 and 5.53 are applied to each of the N nodal points P in turn, in the 


discretised form 


M 
Cor(P)U(P) + Cry (PCP) + Y” J [u(Q)Tax P,Q) + MQITy ,Q)]A5(@ 
m=1 “m 


M 
= | [eee P.O) + fey P,Q] A500) 


m=1 Sm 
(5.54) 

M 

Ce (P)u(P) + Cyy (PCP) + [ [te(Q)T yx (P,Q) + V(Q)Tyy (P,Q) |45(Q) 
m=1 "m 

M 
= f [tx (Q)Uyx (P,Q) + ty (Q)Uyy (P, Q)] dS(Q) 
~ (5.55) 
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The total number of boundary elements is M,m is an element counter, and S,, is the piece of the boundary 
occupied by element number ™. The details of how integration is carried out over an individual element 
depend on the type of element involved. Equations 5.54 and 5.55 represent a set of 2N linear equations, 
where N is the number of nodal points (N = M for constant or linear boundary elements, N = 2M for 
quadratic elements). It is convenient to arrange them in the order of the node numbers, with for each 


node the equation corresponding to the x direction force first, followed by that for the y direction force 
[A][u] = [B][-] (5.56) 


Matrices [A] and [B] are square and contain known constant coefficients, in general all of which will 
be non-zero. The column vectors [u] and [7] contain the nodal point values of displacements and 
tractions, two components for each at each node. At each node in each direction either the displacement 
or traction will be unknown and the other known, or there will be a linear relationship between them. 
The equations can be rearranged, taking all unknown quantities to the left hand side, known quantities 


to the right hand side, giving a set of linear equations in a familiar form 
[A*][x] = [B*]ly] = [4] (5.57) 


where [A*] and [B*] are modified coefficient matrices and [b] is a column vector of known coefficients. 
This set can be solved for the 2N unknowns x at the nodes, meaning that the displacements and tractions 
are known at every nodal point on the boundary. Given this information, stresses at the nodes can also 
be found. 


5.4 Boundary Conditions and Surface Stresses 


The most straightforward type of boundary condition met in practice is where tractions are prescribed. 
In other words, for a particular node the traction components in both of the global co-ordinate directions 
are known. These values contribute to the right hand side of Equations 5.57, the solution of which gives 


the displacement components at the node. 


Traction boundary conditions are most likely to be prescribed in the form of stresses on the boundary, 


in particular the direct stress normal to the boundary and the shear stress along it. 
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Figure 5.7 Stresses and tractions at a point on the boundary 


Figure 5.7 shows the normal stress o,,, and shear stress o,,, acting at a point on the boundary. Note 
the sign convention for shear stresses: a positive shear stress acts in the positive Sdirection, along the 
boundary keeping the domain to the left. The outward unit normal vector 7 at the point concerned is 
inclined at angle y to the x global co-ordinate direction, and has components f, and f, in the x and 
y directions. Because both stresses and both tractions (resultant stresses) act on the same surface, they 
can be resolved like forces 


a 


ty = Onn COSY — Ogn SINY = Onn fy — Osn Ny (5.58) 
ty = Onn SINY + Osn COSY = Onn Ny + Ign Thy (5.59) 


With displacement boundary conditions, for a particular node the displacement components in both 
of the global co-ordinate directions are known. These values must be taken from the left hand to the 
right hand side of Equations 5.57 to contribute to [b]. The solution of the equations gives the traction 


components at the node. 


It should be noted that tractions are often discontinuous at a corner, and indeed at any point where either 
the direction of the outward normal to the boundary or the magnitudes of the stresses, change abruptly. 
Figure 5.8 shows a right-angled corner of a domain with an applied normal stress on one side of the 
corner, and stress free on the other side, a situation that occurs commonly in practice. Label A indicates 
the point of the corner. Consider the two points B, and B, on the two sides of A. At B, the tractions are 
t, = 0,t, = 0, whereas at B, they are ty = 09,t, = 0. These traction values remain unchanged as B, 
and B, approach A. Clearly, there is a discontinuity in traction t, at A, which raises the question as to 


what is meant by the tractions at such a point. 
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Figure 5.8 Tractions at a corner 


In practice, there is no difficulty if stresses are prescribed on both sides of a corner, provided the 
integrals which must be evaluated over both sides of the corner employ the relevant traction values, 
without reference to any (undefined) corner value. The displacement components at the corner, which 
are not discontinuous, are the unknowns to be found. The case of displacements prescribed on both 
sides of a corner is not of great practical interest, corresponding to a form of rigid body displacement. 
At least for uniform displacements the tractions would be zero. If one side of a corner has displacements 
prescribed, and the other stresses, then it is convenient to regard the corner as being part of the prescribed 
displacement side, and treat the unknown tractions as the tractions applicable to this side (the tractions 


on the other side being known). 
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Other types of boundary conditions create greater complexity. For example, a line of symmetry has fixed 
(usually zero) displacement normal to the line, zero traction along it. Such a mixture of displacements 
and tractions can be dealt with, but corners where a line of symmetry meets a part of the boundary with 
either prescribed displacements, prescribed stresses or indeed another line of symmetry, require special 
care. Further possible types of mixed boundary condition include, for example, flexible supports where 


tractions are related to displacements. 


It is important to note that in boundary element methods for elastic stress analysis problems it is not 
possible to apply point force boundary conditions, which would give rise to infinite stresses. This is in 
contrast to other numerical methods such as finite elements. It does, however, reflect physical reality 
in which even concentrated forces on a body are applied over small but finite areas. Similarly, point 
displacement constraints are not permitted if they would require point forces to maintain them. On the 
other hand, such constraints applied to a domain which is already in force equilibrium, simply in order 
to prevent movement of the domain as a rigid body, are permitted. Indeed, they are necessary to prevent 
such a problem giving rise to a singular [A*] matrix in Equations 5.57. This can be a very useful facility, 


and will be demonstrated in the problems considered in Chapter 6. 


Point constraints are straightforward to apply. If, say, a particular node needs to be constrained in the x 
direction, the first of the two equations corresponding to that node is modified. All the coefficients in 
that equation, including that in the right hand side vector, but excluding the coefficient on the diagonal 


of matrix [A*], are set to zero. The diagonal coefficient is set to one. 


The primary results of the analysis are the displacements and tractions at the nodes. Stresses are likely 
to be of much greater interest in practice than tractions. Referring to Figure 5.7, the normal and shear 


stresses can be obtained from the tractions as 
Onn = t, cosy +t, siny =t,fA, + tn, (5.60) 
Osn = —t, Siny + ty cosy = —t, fy + tyn, (5.61) 


The other stress component of interest is the direct stress in the direction along the surface, 95s. This 
cannot be found from the tractions, but can be found from the displacements. If wu, is the displacement 


along the boundary in the S direction, the direct strain there is given by 


ney Ue (5.62) 


from which Oss = E*€s5 +V "Onn (5.63) 
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Given the three stress components, Onn,0,, and o,,,, other stress parameters can also be found, such as 
the principal stresses at the point. Often of greater practical use is a single stress, reflecting the overall 
intensity of the state of stress and hence the likelihood of yielding or failure. The choice of this single or 
equivalent stress therefore depends on the criterion adopted for yielding or failure. At least for ductile 
materials, the most commonly used is the von Mises equivalent stress, which in the present situation 


is defined as 


Oo =V Opn + O02, — Onn Oss + 304, (5.64) 


To test for yield or failure, this equivalent stress is compared with the yield or failure stress for the 


material concerned. 


5.5 Quadratic Boundary Elements 


Experience of potential problems in Chapters 3 and 4 showed that, although constant boundary elements 
are easier to implement, quadratic elements give much more accurate results for the same fineness of 
discretisation. The only exceptions to this were when results were required at points within a solution 
domain close to a boundary, or when the domain was very narrow so that opposite sides of the boundary 
were close together. For elastic stress analysis problems attention is confined here to quadratic boundary 
elements. This is particularly because such problems are often concerned with rapid stress variations in 
local regions such as stress concentrations or cracks. Results at internal points are rarely required, but 


care will need to be taken with long slender domains. 


P 
O 


n(qy 2° 


Figure 5.9 Quadratic element discretisation of a two-dimensional solution domain boundary 
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Figure 5.9 shows a typical arrangement of quadratic line elements on parts of the boundary of a two- 
dimensional solution domain. Each element has three nodes, one at each end and one at its centre. While 
the end nodes are shown as solid circles, those at the centres of elements are shown as open circles. The 
relevant boundary integral equations are Equations 5.54 and 5.55. The total number of nodal points 
is twice the total number of elements, and there are two equations per node, hence 4M equations are 
generated, M being the number of elements. The integrations involved in the evaluation of the equation 


coefficients must in general be carried out numerically. 


Figure 5.10 A typical quadratic line element 
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Figure 5.10 shows a typical curved element, exactly as used for potential problems. Its three nodes are 


numbered in the direction of integration from 1 to 3, these numbers being local to the particular element. 


The second node is at the centre of the element as measured along its curved length. A local intrinsic 


co-ordinate, €, follows the curved shape of the element. It has its origin at the centre node: that is, € = 0 


at node 2, and takes the values -1 and +1 at nodes 1 and 3, respectively. 


The distributions of displacements and tractions are approximated as quadratic functions of position 


along the element in terms of values at the three nodes. For example 
u(&) = Ny (Fu + N2 (Juz + N3(E)Juz = Vea Ne(§) Uc 
Shape functions are as defined in Equations 2.50, 2.51 and 2.52 
M(@) =5éE -1) 
N2(§) =1-& 
N3(6) = 56 + 1) 
Variations of the global co-ordinates are similarly defined 
x(E) = Ny (E)xq + No (E)xz + N3(E)x3 = Vea Ne (E) Xe 
y(E) = Ni @)y + No@y2 + N3@)y3 = Lear Ne (Eve 


The rates of change of the global co-ordinates with respect to ¢ are 


d. 


dx a3 AN (E) dy _ 3 dNc(E) 
é = ees dé Xo> dé — c=1 dé Cc 


and the derivatives of the shape functions are 


aN) _ g 


dN2() dN3(€) _ 1 
A wo) pg, MOA 


d 
Pe dé dé +3 


As in Equation 2.59 the vector normal to the boundary can be found as 


n=n,itnj=SAk=—i-—j 


where n, and ny, are the components in the x and y directions of the vector normal n. 
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(5.65) 


(5.66) 


(5.67) 


(5.68) 


(5.69) 


(5.70) 


(5.71) 


(5.72) 


(5.73) 


In order to integrate round the boundary it is necessary to change from global co-ordinate to local 


intrinsic co-ordinate € within each element, by means of the Jacobian of transformation 


_ ds V@)7+HH_ fax)? , (av\?_ [a 
I@=_= "Re (=) +(2) = |n2 +n2 (5.74) 


which can also be used to define the components of the unit normal 


n ae 


n i ly; 
mae OIG) 


Introducing the transformation from global to intrinsic co-ordinate, the boundary integral equations, 


Equations 5.54 and 5.55, become 


Mgt 
Cox (P)UCP) + Cy PIV) +” | [uCQTx P,Q) + PQMy P, D] JA 
m=1 — 


u +1 
=D J [QQ +4, QU PDO 
m=1 — 


(5.76) 
Cy (P)uCP) + Cyy (PvP) + Ss | “ORs (P,Q) + VQTyy (P.Q)] JES 
Zie 
wv +1 
= 2 [fe @U ys, + & Uy P, ]IE aE es 


Then introducing the parametric representation of Equation 5.65 for the displacement and traction 


distributions 


ee 5 
Cox P)UCP) + Cy PIP) + YY | [ucTec P.O) + veToy (P.QINeCE) JEDAE 
-1 


m=1c=1 


A ae +1 
YY | [ede P.O + {ty} Uo PO] NGO IG 


m=1c= 1 


; (5.78) 
Me Pe ott 
Cyx (P)u(P) + Cyy (P)v(P) + > > i ; cee (P,Q) + UTP Q)|N-(é) J(E)dé 
m=1c=1 — 
Se +1 
2) I, [fe Uyr (PD + {ty} Uyy P,Q] NGG 
m=l1c=1 on 


Download free eBooks at bookboon.com 


Boundary Element Methods for Engineers: 
Part Il: Plane Elastic Problems Boundary Element Analysis of Plane Elastic Problems 


where c is the number (1, 2 or 3) of the node in element number m. The first (traction) kernel functions 
on the left hand sides of the equations and the second (displacement) kernels on the right hand sides 
are given by Equations 5.33, 5.36, 5.38, 5.40 and 5.26 to 5.29. 


5.5.1 Points P and Q not in the same element 


The required integrations are performed using Gaussian quadrature, described in Appendix A. Provided 


point P is not in the element over which integration is required the process is straightforward. 


5.5.2 Points P and Q in the same element 


The case where P is in the relevant element requires some care, because when and are coincident the 
kernel functions are singular. The orders of the singularities are r~ and In(r~*). Provided that P is not 
the cth node of the element, however, the shape function N.(&) goes to zero at P, and the products of 


kernels and shape functions are not singular at P, and may be integrated normally. 


If P is the cth node of the element then the terms at the diagonal of matrix [A], involving both the first 
(traction) kernel functions and the free term constants C,. (P),C, (P),Cy, (P), and C,, (P), are difficult 
to evaluate directly. They can, however, be evaluated indirectly. If there are no tractions applied to the 


boundary of a domain, Equations 5.56 become 


[A]lu] = [0] (5.80) 


22 
a > 
= 
Q 
Qa 
as 
ae 
oO 
ao 
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Possible outcomes of this lack of loading include rigid body displacement of the entire domain in the x 
direction by a unit distance (all x-direction diplacements take the value 1, all y-direction displacements 
take the value 0). Also, rigid body displacement of the entire domain in the y direction by a unit 
distance (all y-direction displacements take the value 1, all x-direction displacements take the value 0). 
This means that the sum of all the coefficients multiplying x-direction displacements in any row of [A] 
should be zero, and, separately, the sum of all the coefficients multiplying y-direction displacements in 
any row of [A] should also be zero. These conditions allow the coefficients of displacements at point P 
(four coefficients in total, two in each of the two equations) to be determined by summations equivalent 
to Equations 2.47 and 2.66. 


Still considering P to be the cth node of the element in Equations 5.78 and 5.79, the second (displacement) 
kernels require special treatment to deal with the In(7r~') singularity, using a modified form of Gaussian 


quadrature. The procedures follow closely those described in Section 2.5.2, but are repeated here for 


convenience. The radial distance from P to a point Q at (x,y) can be arranged in the form 


r(P,Q) = mR) (5.81) 
where R(é) is a known function and n is a modified intrinsic co-ordinate with its origin at P. The form 
of R(¢) depends on which of the three element nodes P is located at. Co-ordinate nis chosen in each case 
to conform to the requirements of Gaussian quadrature involving a logarithmic function (Appendix A). 
P at the first node 
If P is at the first node of the element 

2 

r(P,Q) = [(@-m)’ + -y1)°P (5.82) 

and x and y are defined in terms of intrinsic co-ordinate ¢ by Equations 5.69 and 5.70. So 


x — xy = [Ny (S) — 1]x1 + N2(S)x2 + N3(E)x3 (5.83) 


y— v1 = [Ni (E) — 1]y1 + No (E)y2 + N3(Oy3 (5.84) 
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Intrinsic co-ordinate 7 is chosen to range from 0 at the first node (¢ = —1) to 1 at the third node 
(€ = +1), so that 


n => +1) a =2 (5.85) 

[N, ) - 1] = 5? -&-2) =5 + DE - 2) =n - 2) (5.86) 

N2 (6) = 1-€? = 2n(1 -&) (5.87) 

N3(§) = iG + =e (5.88) 
Therefore 

x — x, = n[(E — 2)x, + 2(1 — E)xq + Ex] (5.89) 

y—yM =n — 2)y1 + 201 — F)y2 + Fys] (5.90) 


R(E) = {[ — 2)x, + 204 — E)xq + Ex3]? + [(E — 2), + 201 — Dy2 + éy,)?}3 (5.91) 


The integrals of the second kernels in Equations 5.78 and 5.79 can be expressed in the form 


+11 +v* eee 


1 
[- Uy (P, Q)N, (E) (EN dE = li 4nE* |(1+v* (Ea Q) 


| + if, 


No (ES) J (ds 


= Kk, [ in(=) ne) [Fan + kk [ow In( a) wel UE + Ki [- APNE — (5,92) 


=) 


+1 (1 +y sf 


+i 
[Yo @on.@s@e = [ 


+1 
AN EIE)E = Ke | FANE IO (5.93) 


tla +y"* y 


+1 
[- Uys PONE IE AE = fT HAN OI@A =Ki [HAN ODI@AES) 


+1 +1 *\2 _ 4y* 
[/ ay tPameersenag = [S| En (I ) + a6] necea seas 
a _,  4mE* (1+v*) \r(P,Q) 


=k f n2\necr@ Elor+ mike [in Gea)ncou@ae + [ RANEIOA (5.95) 


where Ky, soe and K, = Sie (5.96) 
An E* 
are constants. 


Because R(€é) is not zero within the range of integration, all of the integrals except those involving In 


(-) can be evaluated by normal Gaussian quadrature. Those involving the singular logarithmic function 


Ing hea be fxyaluated using the appropriate quadrature formula described in Appendix A. 
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P at the second node 


If P is at the second node of the element 


r(P,Q) = [(—x2)? + (—y2)F2 (5.97) 
X= X_ = Ny (E)xy + [No (SE) — 1]x2 + N3(S)x3 (5.98) 
Y— Yo = Ny (E)y1 + [N2(€) — 1]y2 + N3(E)y3 (5.99) 


For integration purposes, the element needs to be divided into two regions, from the second node to 
the third node (€ = 0 to 1) and from the second node to the first node (€ = Oto —1). Between the 


second and third nodes the intrinsic co-ordinate is chosen as 


m=5 = =A (5.100) 
and Ni (§) = 58 - 1) = 5m - 1) (5.101) 
[No(§) — 1] = -§? = —m§ (5.102) 
N3(€) = 56 +1) =5mE +) (5.103) 
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Therefore 
1 1 
X—-X2 = [5 (é = 1)x, = EX2 + 56 + 1)xs| (5.104) 


y= 92 =I F CDi sy2 eG - 1ys| (5.105) 


NIB 


2 


1 1 1 1 
ne) = {56-2 - 62456 +0] [6 -Di-o2+56 + Dos] : (5.106) 
Between the second and first nodes the intrinsic co-ordinate 7 is chosen as 


ale 
dn2 


ng = -§ (5.107) 


and R(€) is as defined in Equation 5.106. 


The integrals of the second kernels U,,,(P,Q) and U,, (P,Q) in Equations 5.78 and 5.79 are as in 
Equations 5.93 and 5.94. The integrals of the other two can be expressed in the form 


"la +v")* (3=) 
4nE* |(1+4+0*) 


i Uyx (P, QUNe G) J(E)dE = i In (aa D) +h,P | (€) J(@dé 


= ike n(2) neo [Eons + kiko [n(=) ne [Fa 


+1 +1 
tek [ in(g) eae + Kf RANE COIOME (5.108) 
sie ivy (Gav) 1 - 
[, @.on-@s1@ag = | ane ay"(aey) + fy) NCE) 1G 


= ike [in(~) nee [fan + Kite [in( 4) ncev@ [Z| a 


+1 


+1 
bike [ing NEU@aE+K [HAN OIEAE (5.109) 
The third and fourth integrals on the right hand sides of Equations 5.108 and 5.109 can be evaluated by 


normal Gaussian quadrature, while the first and second involve the singular logarithmic function and 


must be evaluated using the appropriate quadrature formula. 
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P at the third node 


If P is at the third node of the element 
2 215 
r(P,Q) = [(«— x3)° + (vy — y3)*]2 
x — X3 = Ny (E)x1 + N2(E)x2 + [N3() — 1] x3 


y —y3 = Ny ()y1 + No(S)y2 + [N3(E) — 1] v3 


(5.110) 


(5.111) 


(5.112) 


Intrinsic co-ordinate 7 is chosen to range from 0 at the third node (¢ = +1) to 1 at the first node 


(& = —1), so that 
n=4a-H |f|= 
M1 (@) = 58 - 1) = -n& 
N2(E) = 1-7 = 2n(1 +8) 
[N3(€) - 1] => 2 + - 2) =F - E+ 2) =-nE +2) 
Therefore 


C= a5 = Hl gag + 20 +S )xe = (E + 2) a5] 


y—y3 = nl-Fy, + 201 + Ey. — E + 2)y3] 


R(®) = {[-Exy + 201 + E)xq — (€ + 2) x3]? + [-Ey. + 21 + E)y2 - 


(5.113) 
(5.114) 
(5.115) 


(5.116) 


(5.117) 


(5.118) 


(E+ 2y,P}7 6.119) 


The integrals of the second kernels U,, (P,Q) and Uy, (P,Q) in Equations 5.78 and 5.79 are as in 
Equations 5.93 and 5.94. The integrals of the other two can be expressed in the form 


Ta +v*)2 B=), 1 
mE* \(1+v*) (aa 


ie Uys (P, Q)Ne(E) J(@)4E = ir 


= ke f'n (2) nce [Elan + Kake [1m (gen) NIUE | 
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)+ nt | Ne) JE) 


+1 (5.120) 


PF No (SI (EE 


= eu Ha +v*)? B-v), it rie : 

[ Uy (PONOIOA = f EE en ( a) +855 -(@) J@)ae 

= Kk { in(2\n G8) nt KK ( in(—)n acak, [ AAn dz (5-121) 
=K, af n(=) neC@ [Fan + af. n (arp) Ne(IG AE + if. £,f,N-(E)I(@)aé 


The second and third integrals on the right hand sides can be evaluated by normal Gaussian quadrature, 
while the first involves the singular logarithmic function and must be evaluated using the appropriate 


quadrature formula. 


5.6 Scaling 


As in the case of two-dimensional potential problems (Section 2.6), it is necessary to scale all distances 
between nodes to be less than unity. This is done by dividing by the maximum dimension of the problem 
(the maximum distance between any pair of nodes). Displacements have the units of length, so these 
are also scaled by dividing by the maximum dimension. Similarly, stresses and tractions are scaled by 
dividing by the value of Young’s modulus. Once the solution to the scaled problem has been obtained, 


the scaling is removed. 


5.7 Solving the Linear Equations 


As for potential problems (Section 2.7), the linear equations arising from the boundary element analysis 
are most appropriately solved by direct methods such as Gaussian elimination. The method is described 


in detail, including an appropriate computer subprogram, in Appendix B. 


5.8 Body Forces 


Many potential problems involve the Poisson form of governing differential equation rather than the 
simpler Laplace form. As a result, a particular integral has to be added, which can be done by modifying 
the boundary conditions for the Laplace problem, as explained in Section 2.8. In elastic stress analysis 
problems the equivalent situation arises if there are significant body forces. The most common of these 
is gravity, but other examples include centrifugal loading. Gravitational effects only become significant 
in physically very large components and structures, such as bridges and dams. In the large majority of 


problems body forces can be neglected, and they are not considered in detail here. 


Problems 


5.1 Under what conditions are the states of plane stress and plane strain indistinguishable? Is this 


of any practical significance? 


5.2 Find the particular forms of the typical displacement and traction kernel functions, 
yx (P,Q) U xy (P1 Q)> Tex (P19) and T,,(p,q) for an incompressible material under plane 


strain conditions. 


Download free eBooks at bookboon.com 


5.3 The typical point on a boundary shown in Figure 5.7 is subject to flexible boundary conditions 
Onn = —KpUn sn > kus 

where u, and u, are the displacements in the n and S directions, and k, and k, are normal 

and shear stiffnesses. Define the boundary conditions for the displacements and tractions in the 


global co-ordinate directions. 


5.4 Assuming that gravity acts in the negative y direction, show that the following distributions of 


stresses and displacements provide a possible particular integral for this body force 


= = * — 
Oyy = PgY, Oxx = V Dyy » Oxy = 0 
pgy* 2 
u=0, ers (1 -—v**) 


where p is the material density, and g is the acceleration due to gravity. 


re If a body whose material has a density of P is rotated at an angular velocity of w about the y 
axis, the centrifugal body force per unit volume generated is pw?x in the x direction. Show that 


the following distributions of stresses and displacements provide a possible particular integral 


for this body force 
Dig d 
pw°x * = 
Oxy = a, Oyy = V' Oxy , Oxy = 0 
2 id 
Wx 
i=- (1 —v*), y=0 
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6 Quadratic Boundary 
Element Program for Plane 
Elastic Problems 


In this chapter a Fortran computer program to implement the quadratic boundary element formulation 
for two-dimensional elastic stress analysis problems developed in Chapter 5 is presented and described 


in detail. It is then used to solve a range of problems to demonstrate the capabilities of the method. 
For readers who prefer to use Matlab, a translation is provided in Appendix E. 


As far as possible the program structure, file names, variable names, and actual coding follow those 
used in the programs for potential problems, particularly the quadratic element program described in 
Chapter 4. As in the case of that program, the principle adopted in programming is to try to make the 
coding straightforward to understand and follow, rather than necessarily the most efficient in terms 


of computation. 
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6.1 Program BEM2EQ 


The program name indicates that it is for boundary element method analysis of two-dimensional elastic 
problems using quadratic elements. Each of the subprogram units which make up the whole are described 


in turn. The Preface explains how the full program can be accessed as a single file. 


As with programs for potential problems, much of the coding is concerned with the definition of the 
arrangement of elements on the boundary (or boundaries) of the solution domain, and the application of 
boundary conditions to them. Again, entering the co-ordinates of the nodes of each and every element, 
followed by the type and magnitude of the relevant boundary condition applied to each is an option, 
but tedious. Instead, each boundary can be divided up into a series of boundary segments, which are 
either straight lines or circular arcs. The number of elements within a segment can then be chosen, and 
varied easily, and from which the program generates all the element geometric data. The elements on 
a segment do not have to be uniform in size, but can be varied in length by a constant ratio between 
successive elements. In the present version of the program, each segment is subject to only one uniform 
boundary condition: defined displacements or defined stresses. The program distributes this condition to 
all the elements involved. Consequently, the ends of segments are conveniently defined as points where 
there is a significant change in either shape (a corner, for example) or boundary condition. Use of this 


facility is demonstrated later in this chapter. 


6.1.1 Main program 


At the beginning of the program is a storage module named SHAREDDATA2EQ which allows stored data 
to be accessed and shared by all those subprograms that require them (by means of a USE statement). 
The dimensioned array sizes in the module allow for up to 250 quadratic boundary elements with 500 
nodes, with up to 10 different boundaries forming the solution domain, and a maximum number of 
point (node) displacement constraints of 20. With two equations per node, this means that up to 1000 
equations can be solved. A dictionary of the variable names used is provided at the beginning of the 


book, and at the beginning of Part II. 


The main program named BEM2EQ is designed mainly to call each of the other subprograms in turn. It 
does, however, also serve to name the files with which the program communicates via OPEN statements. 
File DATA, which is addressed as file number 5 in the program, serves to supply the input data which 
defines the current problem. The main output of results is to file RESULTS, addressed as file number 
6. Element mesh data, on the other hand, are output to file MESHRES (mesh results) and numbered 7. 


MODULE SHAREDDATA2EQ 


! MODULE STORING SHARED DATA. 
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REAL :: XEEND (260) ,YEEND (260) , XNODE (500) , YNODE (500) 
REAL :: XSEND(250),YSEND (250) 
REAL :: UNMX(250,3),UNMY(250,3),ANGSTORE (260) ,MAXL 
REAL :: A(1000,1001),UV(1000) ,U(500) ,V (500) , USEG (250) , VSEG (250) 
REAL :: AROWX (250, 6) ,AROWY (250, 6) ,BROWX (250, 6) , BROWY (250, 6) 
REAL :: ZG(8),WG(8),EGL(8) ,WGL(8) , JACOB, UNGX, UNGY 
REAL :: SIGNNSEG (250) ,SIGSNSEG (250) , UELEM (250) , VELEM (250) 
REAL :: SIGNN(250,3),SIGSS(250,3),SIGSN(250,3),SIGE(250, 3) 
REAL :: TX(500),TY(500), TMX (250, 3) , TMY (250, 3) 
REAL :: FXSEG(250),FYSEG(250) 
REAL :: SF(3,8),SD(3,13),SFL(4,3,8),SDL(4,3,8) 
REAL :: PI,E,NU,ESTORE 
REAL :: AKXX,AKXY,AKYX, AKYY, BKXX, BKXY, BKYX, BKYY 
REAL :: BK2XX,BK2XY,BK2YX, BK2YY 
INTEGER :: NEL,NNP,MAXNEL,MAXNNP,MAXNB, NEEND, NGAUSS 
INTEGER :: NODE (250,3),M1(250) ,M3(250) 
INTEGER :: NBOUND,NSEGTOT, NELB (10) ,NSEGB (10) 
INTEGER :: NBCU, NBCS,NBCM,NBCT, IBCE (250), IBCN (500), ISEGBC (250) 
INTEGER : ISEGEND (250) , ISEGELEM (250) ,MFIRST (250) ,MLAST (250) 
INTEGER :: NBCPC,MAXNPC,NODEPC (20) ,IDIRPC (20) 
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END MODUL 


r 


SHAREDDATA2EQ 


PROGRAM BEM2EQ 


cal 


! PROGRAM FOR SOLVING TWO DIMENSIONAL ELASTIC STRESS ANALYSIS PROBLEMS 


! BY THE BOUNDARY ELEMENT METHOD USING QUADRATIC ELEMENTS. 


USE SHARE 


s 


DATA2EQ 


E="DATA") 


OPEN (5, FI 


OPEN (6, FILE="RESULTS") 


OPEN (7, FILE="MESHRES") 


PI=4.0*ATAN(1.) 


! DEFINE THE MAXIMUM PROBLEM SIZE 


U 


ERMITTED BY THE ARRAY DIMENSIONS. 


MAXNEL=250 
MAXNNP=500 
MAXNB=10 
MAXNPC=20 


| INPUT THE PROBLEM TITLE AND TYP 


Gl 


, ALSO MATERIAL PROPERTIES. 


CALL INTITLE 


! INPUT AND GENERATE THE MESH DATA. 
HO 


CALL M 


Gl 
wn 


! OUTPUT THE 


< 


ESH DATA. 


CALL MSHOUT 


! EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND THEIR DERIVATIV 


= 


Gl 
wn 


! AT THE GAUSS POINTS AND NODES. 


CALL SHAPE 


! INPUT, PROCESS AND OUTPUT THE BOUNDARY CONDITIONS. 


CALL BCS 


! FORM THE COEFFICIENT MATRIX AND APPLY TH 


GI 


BOUNDARY CONDITIONS. 


CALL FRMTRX 
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! SOLVE THE LINEAR EQUATIONS. 


NEQN=2 *NNP 
MAXNEQN=2 *MAXNNP 


MAXNEQNP1=MAXNEQON+1 
CALL ELIMIN (A, UV,NEQN, MAXNEON, MAXNEONP1, IFLAG) 


IF(IFLAG == 1) THEN 


WRITE (6, 61) 


61 FORMAT (/ "MATRIX ILL-CONDITIONING DETECTED IN EQUATION SOLVER") 


STOP 
END IF 


! OUTPUT THE NODAL DISPLACEMENTS, ELEMENT STRESSES AND FORCES ON THE 


! BOUNDARY SEGMENTS. 


CALL OUTPUT 


STOP 


END PROGRAM BEM2EQ 


After defining the maximum numbers of boundary elements (MAXNEL), nodes (MAXNNP), boundaries 
(MAXNB) and point constraints (MAXNPC) permitted by the array dimensions, the main program calls 


the following subprograms in turn: 


INTITLE for the problem title, 

MESHQ to input and create the mesh data, 

MSHOUT to write out the mesh data, 

SHAPE for defining the element shape functions, 

BCS for the boundary conditions, 

FRMTRxX to define the [A] and [B] coefficient matrices, 
ELIMIN to solve the equations, and finally 

OUTPUT to write out the results. 


As in the potential programs, the matrix ill-conditioning warning is most likely to be triggered by trying 
to solve a poorly defined problem, such as one with only prescribed stress boundary conditions, and 
displacement nowhere defined. Computed results from ELIMIN are contained in array UV: typically 


nodal point displacements, but tractions for nodes subject to prescribed displacements. 
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6.1.2 Subprogram INTITLE 


An alphanumeric title for the problem is first read into TITLE. Next a six character message, either 
“STRESS” or “STRAIN” is read into CASE, to define whether the domain is in a state of plane stress or 
plane strain. The default is plane stress: any message other than STRAIN will result in a state of plane 
stress being assumed. Finally, the values of Young’s modulus and Poisson’s ratio are read into E and 
NU. If plane strain is to be assumed, the effective values of these properties are defined according to 


Equations 5.7. 


SUBROUTINE INTITLE 


! SUBPROGRAM TO INPUT PROBLEM TITLE AND WHETHER PLANE 


! STRESS OR PLANE STRAIN. ALSO THE MATERIAL PROPERTIES. 


USE SHAREDDATA2EQ 


CHARACTER(80) :: TITLE 


CHARACTER(6) :: CASE 


| INPUT THE PROBLEM TITLE. 


READ (5, FMT="(A80)") TITLE 


WRITE(6,61) TITLE 
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QuadraticBoundaryElement Program forPlaneElasticProblems 


61 FORMAT ("QUADRATIC BOUNDARY ELEMENT SOLUTION FOR", 
& " TWO DIMENSIONAL ELASTIC PROBLEM" // A) 
! 
! INPUT THE PROBLEM CASE TYPE 
| "STRESS" DEFINES PLANE STRESS 
| "STRAIN" DEFINES PLANE STRAIN 
! THE DEFAULT IS PLANE STRESS. 
READ (5, FMT="(A6)") CASE 
WRITE (6,62) CASE 
62 FORMAT (/"UNDER PLANE ",A," CONDITIONS") 


! INPUT AND OUT 


READ E 


oll 


U 


WRITE E,NU 


63 FORMAT (/"YOUNG'S MODUI 


! MODIFY PROP 


ERTIES IF CASE 


y 
b, 


IF (CASE == "STRAIN") TH 


(1.-NU**2) 


U/ (1.-NU) 


ETURN 


ND SUBROUTINE 


r. 


INTITLI 


Gl 


6.1.3 


PUT YOUNG'S MODULUS AND POISSON'S RATIO. 


E12.4,10X,"POISSON'S RATIO = ",F5.3) 


py 


OF PLAN 


STRAIN. 


Subprogram to input and generate the element mesh data 


Subprogram MESHQ is identical to the subprogram described in Section 4.1.3, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 


SUBROUTIN 


r. 


E SHO 


PROGRAM TO READ IN AND G 


! QUADRATIC ELEM 


G 


NTS. 


Bs 


US ‘DDATA2I 


HA 


Gl 


G 


iN 


BRAT 


THE 


KOM 


ETRIC DATA FOR A MESH OF 


! INPUT THE NUMB 


BER OF S 


KE PARATE 


G 


READ (5, *) 


NBOUND 
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! TEST THE NUMBER OF BOUNDARIES. 


61 


I 


eal 


F(NBOUND < 1 .OR. NBOUND > MAXNB) TH 


gal 
Zz 


WRITE (6,61) NBOUND, MAXNB 


FORMAT (/ "NBOUND =",14,2X,"OUTSIDE 


TU 


ERMITTED RANGE 1 TO",14) 


STOP 
ND IF 


EACH BOUNDARY IN TURN INPUT THE NUMBER OF SEGMENTS. 


ach boundary in turn: DO IBOUND=1,NBOUND 


ELB (IBOUND) =0 
MINB=MMAXB+1 


EAD (5,*) NSEGB (IBOUND) 


SEGTOT=NSEGTOT+NSEGB (IBOUND) 


! TEST THE NUMBER OF SEGMENTS. 


F(NSEGTOT < 1 .OR. NSEGTOT > MAXNEL) THEN 


WRITE (6,62) NSEGTOT,MAXNEL 


62 FORMAT (/ "NSEGTOT =",16,2X,"OUTSIDE PERMITTED RANGE 1 TO",1I6) 
STOP 
END IF 
! 
! INPUT THE CARTESIAN GLOBAL COORDINATES OF THE END POINTS OF THE 
! SEGMENTS. TAKE THE END POINTS CONSECUTIVELY, KEEPING THE DOMAIN 
! TO THE LEFT OF THE DIRECTION OF NUMBERING. 
READ (5,*) (XSEND(ISEND) , YSEND(ISEND) , ISEND=1,NSEGB (IBOUND) ) 


! DEFINE THE FIRST END POINT ON THE CURRENT BOUNDARY. 


XEEND (IEEND) =XSEND (1) 


YEEND (IEEND) =YSEND (1) 


R 
PU 


EACH OF THE SEGMENTS (BETWEEN ENDS 1 AND 2, 2 AND 3, ETC.) 


T TH 


Gl 


RADIUS OF CURVATURE (+VE FOR CONVEX WITH CENTRE OF 


RVATURE INSIDE DOMAIN, -VE FOR CONCAVE), THE NUMBER OF 


EM 


ENTS IN THE SEGMENT, AND THE LENGTH RATIO BETWEEN SUCCESSIVE 
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! ELEM 


=r 


NTS IN THE DIRECTION OF NUMBERING. 


ISEGMAX=NSEGTOT 


ISEGMIN=ISEGMAX-NSEGB (IBOUND) +1 


Each segment in turn: DO ISEG=ISEGMIN, ISEGMAX 


READ(5,*) RSEG, NELSEG, RATSEG 


r 


AND TEST THE NUMBER OF ELEM 


INTS SO FAR. 


ey 
ce 
a 
w 
e) 
Gq 
Zz 
‘Ss 
ll 
Zz 
ml 


ELB (IBOUND) +NELSEG 


F(NEL < 1 .OR. NEL > MAXNEL) THEN 


WRITE (6,63) NEL,MAXNEL 


r 


63 FORMAT (/ "NEL =",16,2X,"OUTSIDE PERMITTED RANGE 1 TO",16) 


! FIRST AND LAST ELEMENTS ON CURRENT SEGMENT. 


MLAST (ISEG) =NEL 


MFIRST (ISEG) =NEL-NELSEG+1 


MMAXB=MMAXB+NELSEG 


ate to mainte- 
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! COORDINATES OF THE FIRST END POINT OF THE 


x 
n 


EGMENT. 


ISEND=ISEG-ISEGMIN+1 


XFIRST=XSEND (ISEND) 


YFIRST=YSEND (ISEND) 


! COORDINATES OF THE LAST END POINT OF THE SEGMENT. 


= 


ISEND=ISEND+1 


IF(ISEG == ISEGMAX) ISEND=1 


XLAST=XSE 


YLAST=YSE 


! GENERATE ELEMENT DATA FOR A STRAIGHT SEGMENT. 


IF(RSEG == 0.) THEN 


Bs 


! DEFINE THE ELEMENT 


IND POINT COORDINATES ON THE SEGMENT. 


Each element in turn: DO M=1,NELSEG 


TEEND=IEEND+1 


ISEGEND (IEEND) =ISEG 


IF(RATSEG == 1.) THEN 


EEND) =XFIRST+ (XLAST-XFIRST) * FLOAT (M) /FLOAT (NELSEG) 


EEND) =YFIRST+ (YLAST-YFIRST) *FLOAT (M) /FLOAT (NELSEG) 


IF(RATSEG /= 1.) THEN 


x 
Zz 
Ss) 

W 


(END) =XFIRST+ (XLAST-XFIRST) * (1.-RATSEG**M) 


& / (1.-RATSEG* *NELSEG) 


YEEND (IEEND) =YFIRST+ (YLAST-YFIRST) * (1.-RATSEG**M) 


& / (1.-RATSEG* *NELSEG) 


END DO Each element in turn 


! DEFINE THE ELEMENT NODES AND COORDINATES. 


ITEEND=IEEND-NELSEG 


Each element in turn: DO IELSEG=1,NELSEG 


TEEND=IEEND+1 


M=MFIRST (ISEG) +IELSEG-1 
T1=2*M-1 


IF(ISEG == ISEGMAX .AND. IELSEG == NELSEG) I3=NODE (MMINB, 1) 
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SEG == ISEGMIN .AND. IELSEG == 1) THEN 


va 
{e) 
Wy UO UD 
ea) 
u 
i 
ll 
K 
Z 
iw) 
a 
Z 
7 
be 


E (13) =XEEND (IEE 


E (13) =YEEND (IEE 


E (12) =0.5* (XNODE (11) +XNODE (T3) ) 


E (12) =0.5* (YNODE (11) +YNODE (T3) ) 


! STORE THE NUMBERS OF THE ADJACENT ELEMENTS. 


zr 


M1 (M) =M-1 


M3 (M) =M+1 


END DO Each element in turn 


END IF 


! GENERATE ELEMENT DATA FOR A SEGMENT IN THE FORM OF A CIRCULAR ARC. 
IF(RSEG /= 0.) THEN 


r 


! LOCATE THE CENTRE OF THE ARC. 


xs 
Ss 
A 
io) 
ll 


(XFIRST+XLAST) /2. 


YMID=(YFIRST+YLAST) /2. 
ALSEG=SOQRT ( (XLAST-XFIRST) **2+ (YLAST-YFIRST) **2) 


> 
U 
bal 


ERP2=RSEG**2- (ALSEG/2.) **2 


IF (ABS (ALPERP2) < 1.E-6*RSEG**2) ALPERP2=0. 


IF(ALPERP2 < -1.E-6*RSEG**2) THEN 


64 FORMAT (/ "DATA ERROR FOR SEGMENT NUMBER",16, 
& / “NOT POSSIBLE TO CREATE A CIRCULAR ARC") 
STOP 
END IF 


ALPERP=SQRT (ALPERP2) 


UVF LX= (XLAST-XFIRST) /ALSEG 


UVFLY= (YLAST-YFIRST) /ALSEG 


FACT=1. 


IF(RSEG < 0.) FACT=-1. 
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XCENT=XMID-ALPERP*UVFLY* FACT 


YCENT=YMID+ALPERP*UVFLX* FACT 


! FIND THE ANGLE SUBTENDED THERE BY THE SEGMENT. 


IF(ALPERP /= 0.) ANGSEG=2.*ATAN (ALSEG*0.5/ALP 


Gl 


RP) 


IF(ALPERP == 0.) ANGSEG=PI 


! DEFINE THE ELEMENT END POINT COORDINATES ON THE SEGMENT. 


ANGFIR=ATAN2 (YFIRST-YCENT, XFIRST-XCENT) 

ANGSTORE (IEEND) =0. 

Each element in turn: DO M=1,NELSEG 

IEEND=IEEND+1 

ISEGEND (IEEND) =ISEG 

IF (RATSEG == .) ANG=ANGSEG* FLOAT (M) /FLOAT (NELSEG) 
IF(RATSEG /= .) ANG=ANGSEG* (1.-RATSEG**M) / (1.-RATSEG* *NELSEG) 
IF(RSEG < 0.) ANG=-ANG 

XEEND (IEEND) =XCENT+ABS (RSEG) *COS (ANGFIR+ANG) 

YEEND (IEEND) =YCENT+ABS (RSEG) *SIN (ANGFIR+ANG) 
ANGSTORE (IEEND) =ANG 

END DO Bach element in turn 


> Apply now 


REDEFINE YOUR FUTURE 
AXA GLOBAL GRADUATE 
PROGRAM 2015 


- © Photononstop 


redefining / standards 


agence cdg 
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! DEFINE THE ELEMENT NODES AND COORDINATES. 


ITEEND=IEEND-NELSEG 


Each element in turn: DO IELSEG=1,NELSEG 


TEEND=IEEND+1 


M=MFIRST (ISEG) +IELSEG-1 


I1=2*M-1 


IF(ISEG == ISEGMAX .AND. IELSEG == NELSEG) I3=NODE (MMINB, 1) 


IF(ISEG == ISEGMIN .AND. IELSEG == 1) THEN 
XNODE (11) =XEEND (IEEND-1) 
YNODE (11) =YEEND (IEEND-1) 
END IF 
XNODE (13) =XEEND (IEEND) 
YNODE (13) =YEEND (IEEND) 
ANG=0 .5* (ANGSTORE (IEEND-1) +ANGSTORE (IEEND) ) 
XNODE (12) =XCENT+ABS (RSEG) *COS (ANGFIR+ANG) 
YNODE (12) =YCENT+ABS (RSEG) *SIN (ANGFIR+ANG) 


! STORE THE NUMBERS OF THE ADJACENT ELEMENTS. 


M1 (M) =M-1 


M3 (M) =M+1 


END DO Each element in turn 


END DO Bach segment in turn 


! ADJACENT ELEMENTS FOR END ELEMENTS OF THE BOUNDARY. 


M1 (MMINB) =MMAXB 
M3 (MMAXB) =MMINB 


END DO Bach boundary in turn 


BEND=IEBEND 


NP=NEL*2 
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! DETERMINE THE MAXIMUM DIMENSION OF THE SOLUTION DOMAIN. 


MAXL=0. 
Each node in turn: DO I=1,NNP 

Each other node in turn: DO J=1,NNP 
DIST=SQRT ( (XNODE (1) -XNODE (J) ) **2+ (YNOD 


EA 


(I) -YNODE (J) ) **2) 


IF(DIST > MAXL) MAXL=DIST 


END DO Each other node in turn 


END DO Each node in turn 


RETURN 


END SUBROUTINE 


< 


E SHO 


6.1.4 Mesh data output subprogram 


Subprogram MSHOUT is identical to the subprogram described in Section 4.1.4, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 


SUBROUTINE MSHOUT 


25 


! SUBPROGRAM TO WRITE OUT THE MESH DATA. 


USE SHAREDDATA2EQ 


! OUTPUT THE NUMB 


Bs 


IRS OF ELEMENTS AND NODES, ALSO THE NODAL 


! COORDINATES. 


WRITE (7,71) NEL,NNP, (I,XNODE (I) ,YNODE(I),I=1,NNP) 


71 FORMAT (/ "GEOMETRIC DATA FOR THE MESH" // 


& 10X, "NUMBER OF ELEMENTS =",16 // 


& 10X, "NUMBER OF NODAL POINTS =",16 // 


& "COORDINATES OF THE NODAL POINTS"// 
& Va Gu I 4 Y ") / 
& 2(16,2E12.4) ) 


! OUTPUT THE ELEMENT NODE NUMBERS. 


ea 


WRITE (7,72) (M, (NODE (M, IN) , IN=1,3) ,M=1,NEL) 


12 FORMAT (/ ‘ELEMENT NODE NUMBERS! // 


& 1X,2(10X, 'M NDL ND2 ND3')/ (2(7X,415))) 


! SCALE THE NODAL POINT COORDINAT! 


Ee 
n 
. 


Each node in turn: DO I=1,NNP 


XNODE (1) =XNODE (I) /MAXL 


Download free eBooks at bookboon.com 


66 


Boundary Element Methods for Engineers: 
Part Il: Plane Elastic Problems QuadraticBoundaryElementProgramforPlaneElasticProblems 


YNODE (1) =YNODE (I) /MAXL 


END DO Each node in turn 


RETURN 


END SUBROUTINE MSHOUT 


6.1.5 Subprogram for defining shape functions 


Subprogram SHAPE is identical to the subprogram described in Section 4.1.6, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 


SUBROUTINE SHAPE 


! SUBPROGRAM TO EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND 


! THEIR DERIVATIVES, AT THE GAUSS POINTS AND NODES. 


USE SHAREDDATA2EQ 


! STORE APPROPRIATE COORDINATES AND WEIGHT FACTORS FOR NORMAL GAUSSIAN 


! QUADRATURE IN ARRAYS XG AND CG, ALSO THOSE FOR LOGARITHMIC FUNCTION 


! INTEGRATION IN XGL AND CGL. 
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NGAUSS=8 
ZG (1) =-0.9602898564 
ZG (2) =-0.7966664774 
ZG (3) =-0.5255324099 
ZG (4) =-0.1834346424 
ZG (5) =-ZG (4) 

ZG (6) =-ZG (3) 
ZG(7)=<26 (2) 

ZG (8) =-ZG (1) 

WG (1) =0.1012285362 
WG (2) =0.2223810344 
WG (3) =0.3137066458 
WG (4) =0. 3626837833 
WG (5) =WG (4) 

WG (6) =WG (3) 

WG (7) =WG (2) 

WG (8) =WG (1) 

EGL (1) =0.013320244 
EGL (2) =0.079750429 
EGL (3) =0.197871029 
EGL (4) =0.354153994 
EGL (5) =0.529458575 
EGL (6) =0.701814530 
EGL (7) =0.849379320 
EGL (8) =0.953326450 
WGL(1)=0.164416605 
WEL)(2)0.237525610 
L (3) =0.226841984 
L (4) =0.175754079 
L(5)=0.112924030 
L(6)=0.057872211 
L(7)=0.020979074 
L (8) =0.003686407 


! NORMAL GAUSSIAN QUADRATURE. 


For each Gauss point in turn: DO IGAUSS=1,NGAUSS 
ZETA=ZG (IGAUSS) 


SF (1, IGAUSS) =0.5*ZETA* (ZETA-1.) 


SF (2, IGAUSS) =1.-ZETA**2 


SF (3, IGAUSS) =0.5*ZETA* (ZETA+1.) 
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SD(1,IGAUSS) =ZETA-0.5 
SD(2, IGAUSS) =-2.*ZETA 


SD (3, IGAUSS) =ZETA+0.5 


END DO For each Gauss point in turn 


=r 


! FOUR CASES OF LOGARITHMIC GAUSSIAN QUADRATURE TO CONSIDER. 


Gl 


| IC=1 - INTEGRATION OVER WHOLE ELEMENT FROM FIRST TO THIRD NODE. 


| IC=2 - INTEGRATION OVER HALF ELEMENT FROM SECOND TO THIRD NODE. 


Gl 


Gl 


| IC=3 - INTEGRATION OVER HALF ELEM 


a 


INT FROM SECOND TO FIRST NODE. 


! IC=4 - INTEGRATION OVER WHOLE ELEMENT FROM THIRD TO FIRST NODE. 


|| 


For each case in turn: DO IC=1,4 


For each Gauss point in turn: DO IGAUSS=1,NGAUSS 


ETA=EGL(IGAUSS) 


IF(IC == 1) ZETA=2.*ETA-1. 
IF(IC == 2) ZETA=ETA 
IF(IC == 3) ZETA=-ETA 
IF(IC == 4) ZETA=1.-2.*ETA 


SFL(IC,1,IGAUSS) =0.5*ZETA* (ZETA-1.) 


SFL(IC,2, IGAUSS) =1.-ZETA**2 


SFL (IC, 3, GAUSS) =0.5*ZETA* (ZETA+1.) 
SDL(IC,1, GAUSS) =ZETA-0.5 


SDL(IC,2, IGAUSS) =-2.*ZETA 


SDL(IC,3, IGAUSS) =ZETA+0.5 


END DO For each Gauss point in turn 


END DO For each case in turn 


! SHAPE FUNCTION DERIVATIVES AT THE NODES, STORED AS THOUGH THEY 


> 
a 
eal 


FOR GAUSS POINTS NUMBERED 11, 12 AND 13. 


Each element node in turn: DO IGAUSS=11,13 


IF(IGAUSS == 11) ZETA=-1. 
IF(IGAUSS == 12) ZETA=0. 
IF(IGAUSS == 13) ZETA=1. 


SD(1,IGAUSS) =ZETA-0.5 
SD(2,IGAUSS) =-2.*ZETA 


SD (3, IGAUSS) =ZETA+0.5 


ND DO Each element node in turn 


eal 


RETURN 


eal 


ND SUBROUTINE SHAPE 
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6.1.6 Subprogram for applying the boundary conditions 


The subprogram BCS serves to apply boundary conditions of either the prescribed displacement or 
prescribed stress types. As already indicated, it is assumed that each segment of elements has a uniform 
boundary condition applied to it, which is an important consideration when defining the segments. Also 
applied are any point displacement constraints. The total numbers of segments subject to the first two 
types of condition, and the number of point constraints, are first read into variables NBCU, NBCS and 
NBCPC, respectively. In the case of prescribed stresses it is only necessary to include those segments 


subject to non-zero stresses. 


Prior to reading in the boundary conditions, the components in the global co-ordinate directions of the 
unit normals at every node are calculated with the aid of subprogram JACOBI. Because the direction of 
the normal to the boundary may be discontinuous at an element end node, particularly at a corner of 
the domain, the calculation is performed for every node of every element, and stored in arrays UNMX 
and UNMY. 
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Storage arrays for normal and shear stresses, SIGNN and SIGSN, and tractions in the global co-ordinate 
directions, TMY and TMY, for every node of every element are set to zero. Array IBCE stores the type 
of boundary condition (1 or 2 for the two types, prescribed displacement or prescribed stress) for each 
element. In the absence of other information, the condition at a node is assumed to be zero stresses. The 


default boundary condition type is therefore set as 2, with the corresponding zero values of stresses for 


the elements being stored in arrays TMX and TMY. 


SUBROUTIN 


Gl 


BCS 


! SUBPROGRAM TO INPUT, 


USE DATA2] 


SHARE 


Gl 


! FIRST FIND THE 


Each element in turn: DO M= 


Each element node in turn: 
IGAUSS=IN+10 

IT=1 

IC=1 

CALL JACOBI (M, IGAUSS, IT, IC) 
UNMX (M, IN) =UNGX 

UNMY (M, IN) =UNGY 

END DO Bach element node in 


END DO Each element in turn 
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61 FORMAT (/ "NBCU =",16,3X,"NBCS =",16,3X,/ "NBCT =",16,3X, 
& "OUTSIDE PERMITTED RANGE 0 TO",I16) 
STOP 
END IF 


IF(NBCPC < 0 .OR. NBCPC > MAXNPC) THEN 


WRITE (6,62) NBCPC,MAXNPC 


62 FORMAT (/ "NBCPC =",14,3X,"OUTSIDE PERMITTED RANGE 0 TO",13) 
STOP 
END IF 
! 
! INITIALISE THE BOUNDARY CONDITION STORAGE ARRAYS. 


Each element in turn: DO M=1,NEL 


IBCE (M) =2 


iru 


Each element node in turn: DO IN=1,3 


SIGNN (M, IN) =0. 
SIGSN (M, IN) =0. 
TMX (M, IN) =0. 
TMY (M, IN) =0. 


END DO Each element node in turn 


END DO Each element in turn 


! INPUT, STORE AND OUTPUT THE PRESCRIBED DISPLACEMENT CONDITIONS. 


IF (NBCU > 0) THEN 


READ (5,*) (ISEGBC (IBCU) , USEG (IBCU) , VSEG(IBCU) , IBCU=1, NBCU) 


WRITE (6, 63) 


63 FORMAT (/ "PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS") 


Bach segment with prescribed displacements: DO IBCU=1,NBCU 


ISEG=ISEGBC (IBCU) 
IF(ISEG < 1 .OR. ISEG > NSEGTOT) TH 


ig 
Zz 


WRITE (6,64) ISEG,NSEGTOT 


64 FORMAT (/ "ISEG = ",16,2X,"OUTSIDE PERMITTED RANGE 1 TO",I6) 


Each element on current segment: DO M=MFIRST(ISEG) ,MLAST (ISEG) 


BCE (M) =1 


T] 
UELEM (M) =USEG (IBCU) 
VELEM (M) =VSEG (IBCU) 


END DO Each element on current segment 


WRITE (6,65) USEG(IBCU) ,VSEG(IBCU) ,MFIRST (ISEG) ,MLAST (ISEG) 
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65 FORMAT (/ "U =",E12.4,3X,"V =",E12.4,3X,"ON ELEMENTS", 14, 3X, 
& "TO",I4) 
END DO Bach segment with prescribed displacements 


END IF 


! INPUT, STORE AND OUTPUT THE PRESCRIBED STRESS CONDITIONS. 


IF(NBCS > 0) THEN 


READ (5,*) (ISEGBC (IBCS) ,SIGNNSEG(IBCS),SIGSNSEG(IBCS), 
& IBCS=1,NBCS) 


WRITE (6, 66) 


66 FORMAT (/ "PRESCRIBED STRESS BOUNDARY CONDITIONS") 


Each segment with prescribed stresses: DO IBCS=1,NBCS 
ISEG=ISEGBC (IBCS) 


F(ISEG < 1 .OR. ISEG > NSEGTOT) THEN 


WRITE (6,64) ISEG,NSEGTOT 


STOP 
END IF 


Each element on current segment: DO M=MFIRST(ISEG) ,MLAST (ISEG) 


IBCE (M) =2 
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! FIND THE TRACTIONS AT TH 


E 


ELEMENT NODES FROM THE PRESCRIBED STRESSES. 


! ALSO STORE THE STRESSES AT THE ELEMENT NODES. 


Each element node in turn: DO IN=1,3 


TMX (M, IN) =SIGNNSEG (IBCS) *UNMX (M, IN) -SIGSNSEG (IBCS) *UNMY (M, IN) 


G 
TMY (M, IN) =SIGNNSEG (IBCS) *UNMY (M, IN) +SIGSNSEG (IBCS) *UNMX (M, IN) 


SIGNN (M, IN) =SIGNNSEG (IBCS) 


SIGSN (M, IN) =SIGSNSEG (IBCS) 


END DO each element node in turn 


END DO Each element on current segment 


WRITE (6,67) SIGNNSEG(IBCS) , SIGSNSEG(IBCS) ,MFIRST(ISEG), 
& MLAST (ISEG) 
67 FORMAT (/ "SIGNN =",F12.4,3X,"SIGSN =",E12.4,3X,"ON ELEMENTS", 
a 14, 3%, "TO", 14) 


END DO Each segment with prescribed stresses 


a 
Z 


D IF 


! ASSEMBLE THE VECTOR OF KNOWN VARIABLES, APPLYING SCALING. 


! FIRST INITIALISE THE NODAL DISPLACEMENTS AND TRACTIONS. 


Each node in turn: DO I=1,NNP 


END DO Each node in turn 


Each element in turn: DO M=1,NEL 


Each element node in turn: DO IN=1,3 


I=NODE (M, IN) 


IBCN (I) =IBC 


Gl 


(M) 


! CHECK THE CONDITION APPLIED TO THE ADJACENT ELEM 


r 


INT FOR AN ELEMENT 


a 
a 


D NODE. IF THE CONDITION IS PRESCRIBED DISPLACEMENTS BUT TH 


ira 


! EXISTING CONDITION AT THE NODE IS NOT, THEN IMPOSE PRESCRIBED 


! DISPLACEMENTS. 


MADJ= 

IF(IN == 1) MADJ=M1 (M) 

IF(IN == 3) MADJ=M3 (M) 

IF(IBCE(MADJ) == 1 .AND. IBCN(I) == 2) IBCN(I)=1 
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! STORE KNOWN VARIABLES. 


E 


IF(IBCN(I) == 1) THEN 
IF(IBCE(M) == 1) THEN 
IF(U(I) == 0.) U(1I)=UELEM (M) /MAXL 
IF(V(I) == 0.) V(I)=VELEM (M) /MAXL 


IF(U(I) /= 0.) U(1I)=0.5* (U(1I) +UELEM (M) /MAXL) 


IF(V(I) /= 0.) V(I)=0.5%* (V(1I) +VELEM (M) /MAXL) 


END IF 
END IF 
IF(IBCE(M) == 2 .AND. IBCN(I) == 2) THEN 
IF(TX(I) == 0.) TX(1I)=TMX(M,IN)/E 
IF(TY(I) == 0.) TY(1I)=TMY(M,IN)/E 
IF(TX(I) /= 0.) TX(I)=0.5* (TX (I) +TMX (M, IN) /E) 
IF(TY(I) /= 0.) TY(I)=0.5* (TY (I) +TMY (M, IN) /E) 
END IF 


ND DO Each element node in turn 


eal 


END DO Each element in turn 


! SCALE YOUNG'S MODULUS. 


ESTORE=E 


rs 


! INPUT, STORE AND OUTPUT NODAL POINT DISPLACEMENT CONSTRAINTS. 


! IDIRPC STORES DIRECTIONS OF THE CONSTRAINTS - 1 FOR X, 2 FOR Y. 


! SUCH CONSTRAINTS SHOULD NOT REQUIRE POINT FORCES TO MAINTAIN THEM. 


IF (NBCU == -AND. NBCPC < 3) THEN 


WRITE (6, 68) 


68 FORMAT (/ "INSUFFICIENT DISPLACEMENT BOUNDARY CONDITIONS", 
& "OR POINT CONSTRAINTS") 
STOP 
END IF 
IF(NBCPC > 0) THEN 


READ (5,*) (NODEPC(IBCPC) ,IDIRPC (IBCPC) , IBCPC=1,NBCPC) 


Each point constraint in turn: DO IBCPC=1,NBCPC 


IF (NODEPC(IBCPC) < 0 .OR. NODEPC(IBCPC) > NNP) THEN 


WRITE (6,69) NODEPC (IBCPC) ,NNP 


69 FORMAT (/ "POINT CONSTRAINT NODE",14," OUTSIDE 


Zz 
Q 


EO TO ",I4) 


STOP 
END IF 


IF(IDIRPC(IBCPC) /= 1 .AND. IDIRPC(IBCPC) /= 2) THEN 
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WRITE(6,70) IDIRPC(IBCPC) , NODEPC (IBCPC) 


70 FORMAT (/ "POINT CONSTRAINT DIRECTION",1I3," FOR NODE",14, 
& " OUTSIDE RANGE 1 TO 2") 
STOP 
END IF 
IF (IDIRPC(IBCPC) == 1) WRITE(6,71) NODEPC(IBCPC) 
71 FORMAT (/ "NODE",14," CONSTRAINED WITH U = 0") 
IF (IDIRPC(IBCPC) == 2) WRITE(6,72) NODEPC (IBCPC) 
72 FORMAT (/ "NODE",14," CONSTRAINED WITH V = 0") 
END DO Bach point constraint in turn 


bel 
Zz 
oO 
= 
ry] 


RETURN 


END SUBROUTINE BCS 


For each segment over which displacements are prescribed, the segment number and values of 
displacements are read into ISEGBC, USEG and VSEG, respectively. Taking each of these segments 
in turn, the elements have their boundary condition type set to 1 (in array IBCE), and the values of 
displacements are stored in arrays UELEM and VELEM. The prescribed displacements and numbers of 


the elements to which they are applied are then written out. 


For each segment over which stresses are prescribed, the segment number and values of normal and shear 
stresses are read into ISEGBC, SIGNNSEG and SIGSNSEG, respectively. Taking each of these segments 
in turn, the elements have their boundary condition type set to 2 (in array IBCE). The values of the 
stresses are both stored as nodal point stresses in arrays SIGNN and SIGNS and converted to tractions 
at the nodes of the elements according to Equations 5.58 and 5.59 and stored in arrays TMX and TMY. 


The prescribed stresses and numbers of the elements to which they are applied are then written out. 


For each node in turn, the known variable values are stored: displacements in arrays U and V, or tractions 
in arrays TX and TY. Displacements are scaled by dividing by MAXL, the maximum domain dimension. 
Tractions are scaled by dividing by Young’s modulus. The type of boundary condition at the node is 
stored in array IBCN, obtained from that for the corresponding element. At the ends of segments, which 
may correspond to physical corners in the problem or may correspond only to changes in boundary 
conditions, this can give rise to an ambiguity in type at a node. At a node at which there is a change from 
displacement to stress boundary condition the former is chosen, which is why in the program there is a 
test for displacement boundary condition on the current element, but stress boundary condition at the 
current node (which can arise from the way the boundary conditions have been applied to the elements). 


The displacement boundary condition type is imposed by setting the relevant value in array IBCN to 1. 
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If at a node linking two segments there is a change, not in boundary condition type but in the values of 
either displacements or stresses, it is appropriate to assign average values to the node and store these as 
the known variables for that node. This is achieved in the program by testing every node to see whether 
both the node and element boundary condition is of the same type. If it is, the known variables are 
taken as the stored values for the element if the known variable values are zero. This will apply to all 
element midside nodes, and to all element end nodes which have not been tested in previous elements. 
If the known variable values are not zero, which can only occur if prescribed values are applied to the 
immediately adjacent element and that conditions at the nodes of the element have already been tested, 
then the average of the present element values and the stored values are re-stored. This works both for 
adjacent elements with different values of displacement or tractions and for the much more common 


case of adjacent elements with the same value prescribed: the common value is re-stored. 


This averaging process is really a detail to give sensible values of the variables in the output of the final 
results. What is really important is that the prescribed displacements and tractions are applied element 
by element in forming the [b] column vector in Equations 5.57, so that each of the adjacent elements 
contributes according to its own prescribed values. This will be explained further in connection with 
subprogram FRMTRX. 


In preparation for using scaled values of the displacement kernels, the value of Young’s modulus is set 
to 1, having stored the actual value in the variable ESTORE. 
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Finally in subprogram BCS the nodal point displacement constraints are read in: the node number into 
NODEPC and the direction of the constraint into IDIRPC. Direction 1 is the x direction and direction 
2 is the y direction. A value of 1 means that the node is constrained not to move in the x direction, 
2 that it cannot move in the y direction. In the absence of any displacement boundary conditions, at 
least 3 (and often only 3) point constraints are required to prevent movement of the domain as a rigid 
body, without at the same time unnaturally restricting deformation of the elastic component under 
consideration. Typically, one node is constrained in both co-ordinate directions, while another which is 
suitably located is constrained to move only in the direction which is radial to the fully constrained node. 
Before the point constraint data are read in, a test is applied to the numbers of displacement boundary 


conditions and point constraints in the subprogram, followed by a warning message where necessary. 


6.1.7 Subprogram to form the coefficient matrix 


Subprogram FRMTRX forms the coefficient matrix and right hand side vector of Equations 5.57, and 
stores them as an extended matrix in array A. The extended matrix is simply [A*] with an extra column 
to contain column vector [b]. Each node in turn is treated as point P, the force point for the fundamental 
solution, then integration is carried out over every element in turn. This integration involves the knowns 
and unknowns at each of the three nodes of each element. Node counters I and J are used for P and 
the current element node, respectively, so that 2*J-1 and 2*] are the numbers of the two rows in [A*] 
corresponding to P, while 2*J-1 and 2*J are the numbers of the two columns in [A*] corresponding to 
the element node. The actual integrations of the kernel function products are carried out in subprogram 
INTKER which is called by FRMTRX. What INTKER does is to compute the integrals of the traction 
and displacement kernel functions required in Equations 5.78 and 5.79. If the knowns and unknowns 
at the th node of the current element are tractions and displacements, then these integral quantities will 
contribute to the [A] and [B] matrices, respectively. Initially they are stored (in subprogram INTKER) 
in the two-dimensional arrays AROWX and AROWY as the coefficients that will form the current pair 
of rows of [A], and BROWX and BROWY as the coefficients that will form the current pair of rows of 
[B], but with the contributions from every node of every element kept distinct. The two rows forming 
the pair for each node are for concentrated force applied at the force point P in the x direction (AROWX 
and BROWX) and y direction (AROWY and BROWY). 


SUBROUTINE FRMTRX 


! SUBPROGRAM TO FORM THE COEFFICIENT MATRIX AND RIGHT HAND SIDE VECTOR, 


! MODIFIED TO SUIT THE BOUNDARY CONDITIONS. 


USE SHAREDDATA2EQ 


! DEFINE THE NUMBER OF COLUMNS IN THE EXTENDED COEFFICIENT MATRIX A. 


JMAX=2*NNP+1 
] 
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QuadraticBoundaryElementProgramforPlaneElasticProblems 


Take each node in turn as P: 


an 


THE 


A(2*I-1,J7)=0. 
A(2*I,J)=0. 


END DO 


AISE 


THE 


B MATRICES. 


Each element in 


Each element 
ROWX (M, 2* IN) =0. 
ROWY (M,2* IN-1) 
ROWY (M,2* IN) =0. 
ROWX (M, 2* IN-1) 
ROWX (M, 2* IN) =0. 


ROWY (M, 2* IN-1) 


ROWY (M,2*IN)=0. 


ND DO 
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UP THE 
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DO J=1,JMAX 
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AIIXX=0. 


AIIXY=0. 


AIIYX=0. 


AIIYY=0. 


Each element in turn: 


Each element node in turn: 


AIIXX=AIIXX-AROWX (M, 2* IN-1) 


Download free eBooks at bookboon.com 


DIAGONAL OF MATRIX A. 


DO M=1,NEL 


79 


DO IN=1,3 
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AIIXY=AIIXY-AROWX (M, 2* IN) 
AIIYX=AIIYX-AROWY (M, 2*IN-1) 
AIIYY=AITYY-AROWY (M, 2*IN) 


END DO Each element node in turn 


END DO Each element in turn 


IF(IBCN(I) == 2) THEN 


A(2*I-1,2*I-1) =AIIXX 
A(2*I-1,2*I) =AIIXY 
A(2*I,2*I-1) =AIIYX 
A(2*I,2*1I) =AILYY 


! INITIALISE THE B*T VECTOR COEFFICIENTS. 


BTX=0. 
BTY=0. 
IF(IBCN(I) == 1) THEN 


BTX=-AIIXX*U (1) -AIIXY*V (I) 
BTY=-AILYX*U (I) ~-AITYY*V (I) 
END IF 
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! APPLY THE BOUNDARY CONDITIONS TO THE CURRENT ROWS OF A AND B, BY 


! CONSIDERING EACH ELEMENT IN TURN. 


Each element in turn: DO M=1,NEL 


Each element node in turn: DO IN=1,3 


J=NODE (M, IN) 


| IF DISPLACEMENTS ARE PRESCRIBED OVER THE ELEMENT, 
B 


! INTERCHANGE THE A AND 


COEFFICIENTS. 


IF(IBCE(M) == 1) THEN 


A(2*I-1,2*J-1) =A (2*I-1,2*J-1) -BROWX (M, 2* IN-1) 
A(2*I-1,2*J)=A(2*I-1, 2*J) -BROWX (M, 2* IN) 
A(2*I,2*J-1)=A(2*1I,2*J-1) -BROWY (M, 2*IN-1) 
A(2*I,2*J) =A(2*1I,2*J) -BROWY (M, 2* IN) 
BTX=BTX-AROWX (M, 2* IN-1) *U (J) -AROWX (M, 2* IN) *V (J) 


BTY=BTY-AROWY (M, 2* IN-1) *U (J) -AROWY (M, 2* IN) *V (J) 
END IF 


| IF STRESSES ARE PRESCRIBED OVER THE ELEMENT, STORE THE A MATRIX 


! COEFFICIENTS, EXCEPT AT A CORNER NODE WHERE PRESCRIB 


r 


iD DISPLACEMENTS 


IF(IBCE(M) == 2) THEN 
BTX=BTX+ (BROWX (M, 2* IN-1) * TMX (M, IN) +BROWX (M, 2* IN) *TMY (M, IN) ) 


& /ESTORE 


BTY=BTY+ (BROWY (M, 2* IN-1) * TMX (M, IN) +BROWY (M, 2* IN) *TMY (M, IN) ) 


& /ESTORE 


IF(IBCN(J) == 2) TH 


Gl 


N 
A(2*I-1,2*J-1) =A (2*I-1, 2*J-1) +AROWX (M, 2* IN-1) 
A(2*I-1,2*J)=A(2*I-1, 2*J) +AROWX (M, 2* IN) 
A(2*I,2*J-1) =A(2*1I,2*J-1) +AROWY (M, 2*IN-1) 
A(2*I,2*J) =A(2*1I,2*J) +AROWY (M, 2* IN) 

END IF 


IF(IBCN(J) == 1) THEN 


BTX=BTX-AROWX (M, 2* IN-1) *U (J) -AROWX (M, 2* IN) *V (J) 
BTY=BTY-AROWY (M, 2* IN-1) *U (J) -AROWY (M, 2* IN) *V (J) 
END IF 
END IF 


END DO Each element node in turn 


END DO Each element in turn 
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! STORE THE B*T COEFFICIENTS AS EXTENSIONS OF MATRIX A. 


A(2*I-1, JMAX) =BTX 


A(2*1I,JMAX) =BTY 


END DO Take each node in turn as P 


! APPLY NODAL POINT DISPLACEMENT CONSTRAINTS. 


IF(NBCPC > 0) THEN 


Each point constraint in turn: DO IBCPC=1,NBCPC 


TROW=2* (NODEPC (IBCPC) -1) +IDIRPC (IBCPC) 
Each column in turn: DO J=1,JMAX 
A(IROW,J)=0. 


END DO Each column in turn 


A(IROW, IROW)=1. 


END DO Bach point constraint in turn 


a 
z 


D IF 


RETURN 


END SUBROUTINE FRMTRX 


At the beginning of the subprogram the current rows of array A and all the coefficients of the AROW and 
BROW arrays are set to zero in anticipation of the assembly process. Then for each boundary element 
in turn subprogram INTKER is called to carry out the integrations. Before any boundary conditions are 
applied, the coefficients of the 2 X 2 submatrix at the diagonal of the [A] matrix, which include the free 
terms, are computed as explained in connection with Equation 5.80. For each of the four coefficients, the 
sum of corresponding values along the current row of matrix [A] must be zero, which allows the value 


at the diagonal to be equated to the sum of the other values with its sign reversed. 


Off-diagonal coefficients for any element end node have two contributions, from the two elements which 
share it, stored in the relevant locations in AROWX and AROWY. The coefficients on the current pair 
of rows of the product of matrix [B*] with the vector of known variables are to be accumulated in BTX 
and BTY. These are set to zero initially in anticipation of the element contributions to be added to them. 
If the boundary condition at point P is of the displacement type, the products of the coefficients at the 
diagonal of [A] with the prescribed displacements there are known and must be moved from the left 
hand side to the right hand side of the equation, and stored in BTX and BTY with their signs changed. 


The boundary conditions applied to each of the elements in turn are now used to determine how the terms 
stored in the AROW and BROW arrays must be added to [A] and [B]. It is perhaps helpful to consider 
the relevant fragments of the equations represented by the current rows of the [A] and [B] matrices, 
the fragments being for the (2*IN-1)" and the (2*IN)" equations for the element, corresponding to the 


IN" node of the element numbered M, as a result of the integration over the element 
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.... }AROWX(M,2*IN-1)x w+ AROWX(M,2*IN)X Yj «0. = 


.... BROWX(M,2*IN-1)x {ty }j +BROWX(M,2*IN)x {ty} + shea (6.1) 


.... + AROWY(M,2*IN-1)x uj+ AROWY(M,2*IN)X Uj «00. = 


.... + BROWY(M,2*IN-1)x {t,}; +BROWY(M,2*IN)x {ty} + peteee (6.2) 


where j is the number of the IN" node of the element (J=>\ NODE(M,IN) in programming terms). These 
parts of a pair of equations are presented in a mixture of program and physical variables in an attempt 


to clarify what is going on. 


Tractions known over the element 


If the tractions over the element are known and the displacements unknown then the above products 
involving BROWX and BROWY, which are known quantities, are added to the right hand side vector 
whose coefficients for the current equation are BIX and BTY. Note that because the traction values 
stored in arrays TMX and TMY have not already been scaled, they are divided by the value of Young’s 
modulus stored in ESTORE. With the displacements at the particular node unknown, on the left hand 
side of the equations the computed AROWX and AROWY values are simply added to the relevant [A] 


matrix coefficients stored in A(I,J). 
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If the displacements at the particular node are known, which implies that the node is an end node at a 
point where the boundary condition changes from displacement type to traction type, then the products 
of the AROWX and AROWY terms with known displacements in Equations 6.1 and 6.2, must be taken 
from the left hand side of the equation to the right hand side and subtracted from BTX and BTY. 


Displacements known over the element 


If the displacements over the element (and hence at all three nodes) are known and the tractions are 


unknown then the terms in Equations 6.1 and 6.2 have to be moved to the opposite sides 
.... -BROWX(M,2*IN-1)x {ty}; - BROWX(M,2*IN)x {ty} + me 7 


.... -AROWX(M,2*IN-1)x uj- AROWX(M,2*IN)X ¥; ..... (6.3) 


.... > BROWY(M,2*IN-1)x {t, }; - BROWY(M,2*IN)x {ty} i+ eee = 


> AROWY(M,2*IN-1)x u- AROWY(M,2*IN)X Yj ..-. — (6.4) 


The products of the AROWX and AROWY terms with the known displacements, which are known 
quantities, are subtracted from the right hand side vector whose coefficients for the current equations 
are BTX and BTY. On the left hand sides of the equations the computed BROWX and BROWY values 


are simply subtracted from the relevant [A] matrix coefficients stored in A(LJ). 


Once assembly of the linear equations for the current point P is complete, the right hand side vector 
coefficients accumulated in BTX and BTY are stored in the last column of extended matrix [A]. The 


assembly process is repeated for all nodes as force point P. 


Finally, any nodal point constraints are applied. As already discussed in Section 5.4, these can only be 
used in positions where they would not give rise to concentrated forces at the points concerned. In other 
words they are used to anchor a domain which is in force equilibrium but would otherwise be free to 
move as rigid body. Point constraints are applied as follows. The relevant equation number to be modified 
is the [2(i — 1) + k]", where i is the node number and k is the direction number of the constraint (1 
for x and 2 for y). All the coefficients in that equation, including the one in the right hand side vector, 
are set to zero. The coefficient on the diagonal of matrix [A] is then set to one. These changes have the 


effect of defining the relevant displacement at the node to be zero. 
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6.1.8 Subprogram for integrating kernel function products over an element 


As indicated in the previous section, the purpose of subprogram INTKER is to compute the integrals 
involving the kernel functions for a particular source point P (node numbered I) and element m, and 
which are required in Equations 5.78 and 5.79. The global co-ordinates of P are first stored in XP and 
YP. Then, for each element node in turn (numbered c in the these equations, IN in the subprogram) the 


node number is stored as J. 


SUBROUTINE INTKER(I,M) 


r. 


! SUBPROGRAM TO INTEGRATE THE KERNEL PRODUCTS FOR A PARTICULAR FORCI 


GI 


! POINT P (INDICATED BY NODE NUMBER I) OVER A PARTICULAR ELEMENT 


! (INDICATED BY ARGUMENT M) BY GAUSSIAN QUADRATURE. 


USE SHAREDDATAZ2] 


ea 
10 


! CONSTANT IN DISPLACEMENT KERNELS MULTIPLYING THE LOGARITHMIC TERM. 


r. 


CL=(1.+NU) * (3.-NU) / (4. *PI*E) 


! COORDINATES OF POINT P. 
XP=XNOD 
YP=YNOD 


! SET UP THE ELEMENT NODE LOOP. 


T. 


Each element node in turn: DO IN=1,3 


J=NODE (M, IN) 


! IF P IS NOT THE CURRENT ELEMENT NODE, USE NORMAL GAUSSIAN QUADRATURE. 


IEA 


IF(I /= J) THEN 


Each Gauss point in turn: DO IGAUSS=1,NGAUSS 


! EVALUATE JACOBIAN AND UNIT NORMAL VECTOR COMPONENTS, ALSO THE KERNELS 


! AT THE PARTICULAR GAUSS POINT FOR NORMAL QUADRATURE OVER THE WHOLE 


|! ELEMENT 
Ti 
ICc=1 
CALL JACOBI (M, IGAUSS, IT, IC) 
CALL KERNEL (XP, YP,M, IGAUSS, UNGX, UNGY) 


! ACCUMULATE THE INTEGRALS. 


SFN=SF (IN, IGAUSS) 
Download free eBooks at bookboon.com 
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AROWX (M, 2* IN-1) =AROWX (M, 2* IN-1) +WG (IGAU 


AROWX (M, 2* IN) =AROWX (M, 


QuadraticBoundaryElementProgramforPlaneElasticProblems 


2* IN) +WG (IGA 


AROWY (M, 2* IN-1) =AROWY (M, 2* IN-1) +WG 


AROWY (M, 2* IN) =AROWY (M, 


2* IN) +WG (IGA 


BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WG 


BROWX (M, 2* IN) =BROWX (M, 


2* IN) +WG (IGA 


BROWY (M, 2* IN-1) =BROWY (M, 2* IN-1) +WG 


BROWY (M, 2* IN) =BROWY (M, 2* IN) +WG (IGA 


END IF 


! IF P IS THE CURRENT ELEMENT NODE, SOME 


Ay 
u 
H 
o 
al 
Zz 


END DO Bach Gauss point in turn 


(IGAU 


(IGAU 
USS) * 


(IGAU 


USS) * 


! P AT THE FIRST NODE OF THE 


EMENT. 


! TERMS INVOLVING NORMAL QUADRATURE. 


Each Gauss point in turn: 


USS) *AKXY*SFN* JACOB 


USS) *AKYY*SFN* JACOB 


BKXY* SFN* JACOB 


LOGARITHMIC QUADRATUR 


SS) *AKXX*SFN* JACOB 


SS) *AKYX*SFN* JACOB 


SS) *BKXX*SFN* JACOB 


SS) *BKYX*SFN* JACOB 


BKYY* SFN* JACOB 


Gl 


DO IGAUSS=1,NGAUSS 
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IT=1 
IC=1 
CALL JACOBI (M, IGAUSS,IT,IC) 
CALL KERN2 (M, IN, IGAUSS) 


SFN=SF (IN, IGAUSS) 
BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WG (IGAUSS) *BK2XX* SFN* JACOB 
ROWX (M, 2* IN) =BROWX (M, 2* IN) +WG (IGAUSS) *BK2XY*SFN* JACOB 


B 
BROWY (M, 2* IN-1) =BROWY (M, 2* IN-1) +WG (IGAUSS) *BK2YX*SFN* JACOB 
B 


ROWY (M, 2* IN) =BROWY (M, 2* IN) +WG (IGAUSS) *BK2YY*SFN* JACOB 


eal 


ND DO Each Gauss point in turn 


! TERMS INVOLVING LOGARITHMIC QUADRATURE. 


Ey 


Each Gauss point in turn: DO IGAUSS=1,NGAUSS 
IT=2 

IC=1 

CALL JACOBI (M, IGAUSS,IT,IC) 
SFN=SFL (IC, IN, IGAUSS) 


DZDE=2. 


1B 


BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WGL (IGAUSS) *CL*SFN* JACOB* DZD 


eA 


BROWY (M, 2* IN) =BROWY (M, 2* IN) +WGL (IGAUSS) *CL*SFN* JACOB* DZD! 


END DO Each Gauss point in turn 


ea 
Pan 


D IF 


! P AT THE SECOND NODE OF THE ELEMENT. 


1 


! TERMS INVOLVING NORMAL QUADRATURE. 


Each Gauss point in turn: DO IGAUSS=1,NGAUSS 


IT=1 

Ic=1 

CALL JACOBI (M, IGAUSS, IT,IC) 
CALL KERN2 (M, IN, IGAUSS) 


SFN=SF (IN, IGAUSS) 
BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WG (IGAUSS) *BK2XX* SFN* JACOB 
ROWX (M, 2* IN) =BROWX (M, 2* IN) +WG (IGAUSS) *BK2XY*SFN* JACOB 


B 
BROWY (M, 2* IN-1) =BROWY (M, 2* IN-1) +WG (IGAUSS) *BK2YX* SFN* JACOB 
B 


ROWY (M, 2* IN) =BROWY (M, 2* IN) +WG (IGAUSS) *BK2YY*SFN* JACOB 


eal 


ND DO Each Gauss point in turn 
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! TERMS INVOLVING LOGARITHMIC QUADRATURE. 


Each Gauss point in turn: DO IGAUSS=1,NGAUSS 
IT=2 

IC=2 

CALL JACOBI (M, IGAUSS,IT,IC) 
SFN=SFL (IC, IN, IGAUSS) 

DZDE=1. 


GI 


BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WGL (IGAUSS) *CL*SFN* JACOB* DZD 


BROWY (M, 2* IN) =BROWY (M, 2* IN) +WGL (IGAUSS) *CL*SFN* JACOB* DZDE 
IC=3 

CALL JACOBI (M, IGAUSS,IT,IC) 

SFN=SFL (IC, IN, IGAUSS) 


DZDE=1. 


iB 


BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WGL (IGAUSS) *CL*SFN* JACOB* DZD 


1B 


BROWY (M, 2* IN) =BROWY (M, 2* IN) +WGL (IGAUSS) *CL*SFN* JACOB* DZD! 


END DO Each Gauss point in turn 


ea 
Pa 


D IF 


! P AT THE THIRD NODE OF THE ELEMENT. 


IF(IN == 3) THEN 


Gl 


! TERMS INVOLVING NORMAL QUADRATURE. 


Each Gauss point in turn: DO IGAUSS=1,NGAUSS 


IT=1 

ICc=1 

CALL JACOBI (M, IGAUSS,IT,IC) 
CALL KERN2 (M, IN, IGAUSS) 


SFN=SF (IN, IGAUSS) 

BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WG (IGAUSS) *BK2XX* SFN* JACOB 
BROWX (M, 2* IN) =BROWX (M, 2* IN) +WG (IGAUSS) *BK2XY*SFN* JACOB 
BROWY (M, 2* IN-1) =BROWY (M, 2* IN-1) +WG (IGAUSS) *BK2YX* SFN* JACOB 


BROWY (M, 2* IN) =BROWY (M, 2* IN) +WG (IGAUSS) *BK2YY*SFN* JACOB 


END DO Each Gauss point in turn 


Gl 


! TERMS INVOLVING LOGARITHMIC QUADRATURE. 


Each Gauss point in turn: DO IGAUSS=1,NGAUSS 
IT=2 

ICc=4 

CALL JACOBI (M, IGAUSS,IT,IC) 
SFN=SFL (IC, IN, IGAUSS) 
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DZDE=2. 
BROWX (M, 2* IN-1) =BROWX (M, 2* IN-1) +WGL (IGAUSS) *CL*SFN* JACOB* DZD 


ea 


Es 


BROWY (M, 2* IN) =BROWY (M, 2* IN) +WGL (IGAUSS) *CL*SFN* JACOB* DZD! 


END DO Bach Gauss point in turn 
END IF 
END IF 
END DO Bach element node in turn 


RETURN 


END SUBROUTINE INTKER 
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P not at the current node of the element containing Q 


If P and Q are not in the same element, or P is not at the current node of the element containing Q 
(that is, I#J) then normal Gaussian quadrature (Appendix A) over the whole element is carried out to 
evaluate the integrals. The required shape functions at the Gauss points have already been computed in 
subprogram SHAPE and stored in array SE. The Jacobians at the Gauss points are found in subprogram 
JACOBI, and the kernel functions in subprogram KERNEL. The last of these subprograms stores the 
kernel functions in variables AKXX, AKXY, AKYX, AKYY and BKXX, BKXY, BKYX, BKYY, respectively. 
These correspond to T,,, Ty, Tyx, 
the products of Gaussian weighting factor (WG), kernel functions, shape function (SFN) and Jacobian 
(JACOB) are then added to either arrays AROWX, AROWY, BROWX or BROWY to complete the 


numerical integration process. 


Tyy and Uy, Uxy, Uyx»Uy, in the equations. For each Gauss point 


P at the current node of the element containing Q 


If P and Q are in the same element and P is at the current node of the element (that is, [=J) then, 
because of the singular nature of the second kernel function, some logarithmic Gaussian quadrature is 
required, as described in Section 5.5.2. The first kernel functions do not have to be evaluated because the 
relevant coefficients at the diagonal of the [A] matrix have already been found indirectly by summing 


the off-diagonal terms. 


If P is the first node of the element the integration process follows Equations 5.82 to 5.96. The integrals 
of the second kernel function are split into singular and non-singular parts. To the latter are applied 
Gaussian quadrature in the usual way, the non-singular parts of the second kernel functions being defined 
in subprogram KERN2 and returned in variables BK2XX, BK2XY, BK2YX and BK2YY. The singular 
parts are integrated using logarithmic Gaussian quadrature (Appendix A) applied to the whole element, 
with the origin of the modified intrinsic co-ordinate y being at the first node. The derivative of § with 
respect to 7, |d&/dy], is stored in variable DZDE. 


If P is the second node of the element the integration process follows Equations 5.97 to 5.109. The 
integrals of the second kernel functions are again split into singular and non-singular parts. To the latter 
are applied Gaussian quadrature in the usual way, the non-singular parts of the second kernel function 
being defined in subprogram KERN2 and returned in variables BK2XX, BK2XY, BK2YX and BK2YY. 
The singular parts are integrated using logarithmic Gaussian quadrature applied to the two halves of the 


element separately, in each case with the origin of the modified intrinsic co-ordinate7 at the second node. 


If P is the third node of the element the integration process follows Equations 5.110 to 5.121, and is 


very similar to the case of P at the first node of the element. 


6.1.9 Subprogram for computing the Jacobian at a point on an element 


Subprogram JACOBI is identical to the subprogram described in Section 4.1.10, with the exception that 
the USE statement is for module SHAREDDATA2EQ rather than SHAREDDATA2PQ. 
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SUBROUTINE JACOBI (M, IGAUSS,IT,IC) 


! SUBPROGRAM TO EVALUATE THE JACOBIAN AND THE COMPONENTS OF THE UNIT 


r. 


Tr. 


! NORMAL VECTOR AT A PARTICULAR GAUSS POINT. 


E 


! M INDICATES THE ELEMENT NUMBER. 


! IGAUSS INDICATES THE 


zr 


GAUSS POINT NUMB 


ea 
w 


Gl 


! IT INDICATES THE TYPE OF QUADRATURE. 


| IT=1 - NORMAL GAUSSIAN QUADRATURE. 


| IT=2 - LOGARITHMIC GAUSSIAN QUADRATURE. 


1B 


! IC INDICATES THE CASE NUMBER FOR LOGARITHMIC GAUSSIAN QUADRATUR 


! AS DEFINED IN SUBROUTINE SHAPE. 


USE SHAREDDATAZ2] 


r 


eal 
10 


! CALCULATE THE DERIVATIVES OF THE GLOBAL COORDINATES WITH RESPECT TO 


! THE LOCAL INTRINSIC COORDINATE. 


DX=0. 
DY=0. 


Each element node in turn: DO IN=1,3 
IF(IT == 1) SFND=SD(IN, IGAUSS) 

IF(IT == 2) SFND=SDL(IC, IN, IGAUSS) 
J=NODE (M, IN) 


DX=DX+SFND*XNOD 


GI 
en 

Cy 
— 


DY=DY+SFND*YNOD 


GI 
a 
Cy 
= 


END DO Each element node in turn 


! COMPONENTS OF THE LOCAL OUTWARD NORMAL VECTOR AT THE GAUSS POINT. 


r 


UNGX=DY 


UNGY=-DX 


!  JACOBIAN OF THE COORDINATE TRANSFORMATION. 


JACOB=SORT (UNGX* *2+UNGY**2) 


! SCALE THE VECTOR COMPONENTS TO GIVE THE UNIT NORMAL VECTOR. 


UNGX=UNGX/ JACOB 
UNGY=UNGY / JACOB 


END SUBROUTINE JACOBI 
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By products of the calculation for a Jacobian (JACOB) are the components of the outward normal vector 
to the boundary at the chosen Gauss point. These are converted to components of the unit normal and 
stored in variables UNGX and UNGY. 


6.1.10 Subprogram for computing the kernel functions when P is not at the current node of the 
element containing Q 


Subprogram KERNEL computes the kernel functions at a particular Gauss point Q. The co-ordinates 
XQ and YQ of the point are first found using Equations 5.69 and 5.70, followed by the components 1, 
and 7, (RX and RY) of the radius vector of length r (R) joining P and Q. The rate of change of radius 
r with the normal n to the boundary (DRDN) is found as the cosine of the angle between them, or the 
scalar product of the two unit vectors. The first kernel functions are computed using Equations 5.33, 
5.36, 5.38, 5.40, and the second kernels from Equations 5.26 to 5.29. 


x 
Ds 
Zz 


SUBROUTINE EL (XP, YP,M, IGAUSS, UNX, UNY) 


! SUBPROGRAM TO EVALUATE THE KERNELS WHEN P IS NOT THE CURRENT 


! ELEMENT NODE. 


! XP, YP INDICATE THE GLOBAL COORDINATES OF POINT P. 


! M INDICATES THE ELEMENT NUMBER. 


29 


! IGAUSS INDICATES THE NUMBER OF THE GAUSS POINT, QO. 
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USE SHAREDDATAZ2] 


eal 
10 


! COORDINATES OF GAUSS POINT Q. 
XQ=0. 
YO=0. 


Each element node in turn: DO IN=1,3 
SFN=SF (IN, IGAUSS) 
J=NODE (M, IN) 


XQ=XQO+SFN*XNOD 


GI 
a 

Cy 
— 


YOQ=YOQ+SFN* YNOD 


GI 
a 

qy 
= 


END DO Each element node in turn 


I 


! COMPONENTS AND MAGNITUDE OF THE RADIUS VECTOR FROM P TO Q. 


RX=XQ-XP 
RY=YQ-YP 
RSQ=RX**2+RY**2 
R=SOQRT (RSQ) 


RX=RX/R 
RY=RY/R 


r 


Bs 


zr 


! RATE OF CHANGE OF R WITH THE NORMAL TO THE BOUNDARY AT Q. 


DRDN=RX* UNX+RY*UNY 


! PARAMETERS IN TH 


zr 


KERNEL FUNCTIONS. 


Cl=-1./(4.*PI*R) 
C2=1.-NU 
C3=2.* (1.+NU) 

C4=(1.+NU) **2/ (4. *PI*! 


Gl 


) 
C5=(3.-NU) *ALOG(1./R) / (1.+NU) 


! EVALUATE THE KERNELS. 
AKXX=C1* (C2+C3*RX*RX) *DRDN 

ERM1=C3*RX*RY*DRDN 

ERM2=C2* (RX*UNY-RY*UNX) 


KXY=C1* (TERM1-TERM2) 


KYY=C1* (C2+C3*RY*RY) *DRDN 


KXX=C4* (C5+RX*RX) 


T 
T 
A 
AKYX=C1* (TERM1+TERM2) 
A 
B 
B 


KXY=C4*RX*RY 
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BKYX=BKXY 


BKYY=C4* (C5+RY*RY) 


RETURN 


eal 


ND SUBROUTINE KERNEL 


6.1.11. Subprogram for computing the second kernel function when P is at the current node of 
the element containing Q 


Subprogram KERN2 computes, for a particular Gauss point Q, the non-singular components of the 
second kernel functions when force point is at the current node of the element containing Q. In other 
words, all but the In (-); In (—), and In (<) terms in Equations 5.92 to 5.96, 5.108 and 5.109, 5.120 and 
5.121. These are stored in BK2XX, BK2XY, BK2YX and BK2YY in the subprogram. The value of function 
R(&) (REN in the subprogram) is defined by either Equation 5.91, 5.106 or 5.119, according to whether 
P is at the first, second or third node of the element. The relevant value of the intrinsic co-ordinate € 


(variable ZETA in the program) is obtained from array ZG for normal Gaussian quadrature. 


SUBROUTINE 


~ 


-RN2 (M, IN, IGAUSS) 


UBPROGRAM TO EVALUATE THE NON-SINGULAR LOGARITHMIC TERM IN TH 


Gl 


ECOND KERNEL WHEN P IS THE CURRENT ELEMENT NODE. 


INDICATES THE ELEMENT NUMBER. 


N INDICATES THE NUMBER OF THE ELEMENT NODE FORMING P. 


T. 


H HH ZF WN Ww 


GAUSS INDICATES THE GAUSS POINT NUMBER. 


USE SHAREDDATA2EQ 


! COORDINATES OF GAUSS POINT Q. 


XQ=0. 
YO=0. 


Each element node in turn: DO IC=1,3 
SFN=SF (IC, IGAUSS) 
J=NODE (M, IC) 


XQ=XQ+SFN*XNODE (J) 


YOQ=YQ+SFN*YNODE (J) 


El 


ND DO Bach element node in turn 


! ELEMENT NODE NUMBERS. 


I1=NODI 


EI 


(M, 1) 
I2=NOD! 


Gl 


(M, 2) 
(M, 3) 


I3=NODI 


El 
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Gl 


! EVALUATE THE INTRINSIC COORDINATE. 


ZETA=ZG (IGAUSS) 


! P AT THE FIRST NODE OF THE ELEMENT. 


IF(IN == 1) THEN 


GI 


XCOMP= (ZETA-2.) *XNODE (I1)+2.* (1.-ZETA) *XNODE (12) +ZETA* XNODE (13) 


Gl 


YCOMP= (ZETA-2.) *YNODE (I1)+2.* (1.-ZETA) *YNODE (12) +ZETA* YNODE (13) 
RFN=SOQRT (XCOMP* *2+YCOMP** 2) 
I=I1 


END IF 


! P AT THE SECOND NODE OF THE ELEMENT. 


IF(IN == 2) THEN 
XCOMP=-0.5* (ZETA-1.) *XNODE (11) +ZETA*XNODE (12) 


& -0.5* (ZETA+1.) *XNO 


D 
D 
YCOMP=-0.5* (ZETA-1.) *YNODE (11) +ZETA*YNODE (12) 
& -0.5* (ZETA+1.) *YNOD 


RFN=SOQRT (XCOMP* *2+YCOMP**2) 
I=1I2 
END IF 


‘aM UID. EL0ZO 


z= 
= 
= 
Ga 
a 
G 
2 
a 
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! P AT THE THIRD NODE OF THE ELEMENT. 


XCOMP=- ZETA* XNODE (11) +2.* (1.+ZETA) *XNODE (12) - (ZETA+2.) *XNODE (13) 


YCOMP=-ZETA*YNODE (11) +2.* (1.+ZETA) *YNODE (12) - (ZETA+2.) *YNODE (13) 


RFN=SQORT (XCOMP* *2+YCOMP**2) 
I=I3 
END IF 


! COMPONENTS AND MAGNITUDE OF THE RADIUS VECTOR FROM P TO Q. 


XP=XNODE (1) 
YP=YNODE (1) 
RX=XQ-XP 
RY=YQ-YP 
RSQ=RX**2+RY**2 


R=SORT (RSQ) 


RX=RX/R 
RY=RY/R 


! PARAMETERS IN THE KERNEL FUNCTIONS. 
C4=(1.+NU) **2/ (4. *PI*E) 


C5=(3.-NU) *ALOG(1./REFN) / (1.+NU) 


! EVALUATE THE KERNELS. 


BK2XX=C4* (C5+RX*RX) 


K2XY=C4*RX*RY 


B 
BK2YX=BK2XY 
BK2YY=C4* (C5+RY*RY) 


RETURN 


ea 


ND SUBROUTINE KERN2 


6.1.12 Results output subprogram 


Subprogram OUTPUT organises and writes out the primary results of the computation. Firstly, 


displacements and tractions are stored in their proper arrays. 


SUBROUTINE OUTPUT 


E 


! SUBPROGRAM TO WRITE OUT THE NODAL POINT VALUES OF DISPLACEMENTS 


! AND ELEMENT STRESSES AND COMPUTE THE FORCES ON THE BOUNDARY SEGMENTS. 
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USE SHAREDDATA2EQ 


! ARRANGE FOR U, V AND TX, TY TO CONTAIN THE NODAL DISPLACEMENTS 


! AND TRACTIONS, RESPECTIVELY. 


Each node in turn: DO I=1,NNP 


IF(IBCN(I) == 1) THEN 
TX (I) =UV (2*I-1) 
TY (I) =UV (2*T) 


a 
a 


D IF 


V (I) =UV (2*T) 


END DO Each node in turn 


! HEADING FOR OUTPUT OF NODAL DISPLACEMENTS AND TRACTIONS. 


WRITE (6, 61) 


61 FORMAT (/ "NODAL POINT DISPLACEMENTS AND TRACTIONS" // 


& I U V TX EY) 


E=RBSTORE 


Each node in turn: DO I=1,NNP 


! REMOVE THE SCALING. 


U (I) =U (I) *MAXL 


V (I) =V (1) *MAXL 


TX (I) =TX (I) *ESTORE 


TY (I) =TY (I) *ESTORE 


r 


! OUTPUT THE NODAL VALUES OF DISPLACEMENTS AND TRACTIONS. 


WRITE (6,62) I,U(I),V(I),TX(I),TY(T) 
62 FORMAT (14, 4E15.6) 


END DO Each node in turn 


! HEADING FOR OUTPUT OF ELEMENT STRESSES. 
WRITE (6, 63) 


63 FORMAT (/ "STRESSES AT THE NODES OF THE ELEMENTS" // 


& M IN SIGNN SIGSS SIGSN", 


& " SIGE") 
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! NORMAL AND SHEAR STRESSES AT THE NODES OF THE ELEMENTS. 


Each element in turn: DO M=1,NEL 


Each element node in turn: DO IN=1,3 


IF(IBCE(M) /= 2) THEN 


J=NODE (M, IN) 

TMX (M, IN) =TX (J) 

TMY (M, IN) =TY (J) 

SIGNN (M, IN) =TMX (M, IN) *UNMX (M, IN) + TMY (M, IN) *UNMY (M, IN) 


SIGSN (M, IN) =-TMX (M, IN) *UNMY (M, IN) +TMY (M, IN) *UNMX (M, IN) 


END IF 
! 
! DIRECT STRESS ALONG THE BOUNDARY SURFACE. 
IC1=NODE (M, 1) 
IC2=NODE (M, 2) 
IC3=NODE (M, 3) 


IGAUSS=IN+10 
DUDZ=SD (1, IGAUSS) *U(IC1)+SD(2, IGAUSS) *U (IC2)+SD(3, IGAUSS) *U(IC3) 


DVDZ=SD (1, IGAUSS) *V(IC1) +SD(2, IGAUSS) *V (IC2)+SD(3, IGAUSS) *V(IC3) 
IT=1 
Ic=1 
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CALL 


SIGSS (M, IN) =E 


VON MIS 
SIGI 


& 


OUTPUT TH 


WRI 


& 


64 


ACC 


COM 


APP 


* 


G 


ES 


KQUIVALI 


ENT STR 
E (M, IN) =SORT (SIGNN 


QuadraticBoundaryElementProgramforPlaneElasticProblems 


JACOBI (M,10+IN, IT, IC) 


ESS. 
(M, IN) 


ESS= (-DUDZ*UNMY (M, IN) +DVDZ*UNMX (M, IN) ) / (JACOB*MAXL) 
ESS+NU*SIGNN (M, IN) 


**2+SIGSS (M, IN) **2 


-SIGNN (M, IN) *SIGSS (M, IN) +3. *SIGSN (M, IN) **2) 


r 


STRE 


ae 


E (6, 64) 


FORMAT (214, 4E15.6) 


END 
END 


EB 
ah 


PUTE 


DO 
DO 


E THE FORC 


FXS 


EG (ISEG)=0. 


FYS 


ND 


ac 


EG (ISEG)=0. 


DO 


h elemen 


a 


G=ISEGEL 


EM (M 


G 


FXE 


1EM=0 . 


G 


FYE 


1EM=0 . 


Each element in 


ES ON TH 


ach segment in turn: 


in turn: 


) 


Each Gauss point in turn: 


SFN=SF (IN, IGAUSS) 


IT= 
[C= 


1 
1 


CAL 


}LEM=FXEL 


G 


JEM=FYEL 


DO 


DO 


HE 


G 


NOD 


EL 


G 


ES OF THE 


EM 


ENTS. 


N(M, IN 


E (M, IN) 


turn 


GI 


BOUN 
DO IS 


Each segment in turn 


DO M= 


LY GAUSSIAN QUADRATURE. 


1 | 


Each element node in turn: 


D 


EM+WG (IGAUSS) *SFN*J 
EM+WG (IGAUSS) *SFN*J 
Each Gauss point in tu 


Each element node in t 


), SIGSS (M, IN) ,SIGSN(M, IN), 


Each element node in turn 


DARY S 


EGM 


ENTS. 


EG=1,NSEGTOT 


Gl 


1,N 


DO IN=1,3 
O IGAUSS=1,NGAUSS 


, JACOBI (M, IGAUSS,IT,IC) 


ACO! 


B* TMX (M, IN) *MAX] 


ACOB* TMY (M, IN) *MAX] 
cn 


ETT 


ENT. 
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FYSEG (ISEG) =FYSEG (ISEG) +FYELEM 


END DO Each element in turn 


! OUTPUT THE SEGMENT FORCE RESULTS. 


WRITE(6,65) (ISEG, FXSEG(ISEG) , FYSEG (ISEG) , ISEG=1,NSEGTOT) 


65 FORMAT (/ "FORCES ACTING ON THE BOUNDARY SEGMENTS" 
& // “SEGMENT FX FY" / (15,2E14.5)) 


RETURN 


END SUBROUTINE OUTPUT 


If the boundary condition at a particular point is of the first type (prescribed displacements) the solution 
array UV contains the computed tractions, which are stored in arrays TX and TY. If the boundary 
condition is of the second type then displacements were the unknowns and are now stored in arrays 
U and V. Then all displacements and tractions have the scaling removed, and nodal point values of 


displacements and tractions are written out. 


Stresses at each of the three nodes of each element are then computed, unless they were given as input 
data. In many cases this will provide duplicate information for end nodes of adjacent elements, but not 
if such nodes are at corners of the domain where tractions are discontinuous. Normal and shear stresses 
are first found using Equations 5.60 and 5.61, followed by the direct stresses along the boundary from 
Equations 5.62 and 5.63. The von Mises equivalent stress is found from Equation 5.64, and all four 


stresses are written out. 


In many problems of practical engineering interest, it is useful to have the global co-ordinate components 
of the forces applied to each of the boundary segments. These can found by integrating the tractions 
element by element (forces stored in FKELEM and FYELEM), summing to give total segment forces, 
FXSEG and FYSEG, which are then written out. 


6.1.13 Subprogram for solving the linear equations 


Both the Gaussian elimination algorithm and subprogram ELIMIN for solving the linear equations are 
described in Appendix B. The matrix [A*] extended to include the right hand side vector [b](in Equation 


5.57) is entered in array A, and the solutions returned in array X (UV in the calling main program). 


6.2 Some Test Problems for BEM2EQ 


Before attempting to solve real problems for which the solutions are not known, it is essential to test 
any computer program as extensively as possible on problems for which exact analytical solutions are 
available for comparison. Test problems serve not only to verify that the program is working correctly, 
but also to examine the level of discretisation (the fineness of the boundary element mesh in this case) 


necessary to obtain accurate results. 
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6.2.1 Simple tension 


One of the most basic test problems available is simple tension. Figure 6.1 shows a thin strip of material 
which has dimensions of 4 units by 2 units. Both the left and right hand edges of the strip (AD and 
BC) are subject to uniform normal tensile stresses of magnitude 3 units. The top and bottom edges (DC 
and AB) are free of applied stresses. Young’s modulus for the material of the strip is 2000 units, and its 
Poisson's ratio is 0.3. Under plane stress conditions, with a uniaxial stress of o,, = 3 the direct strain in 
the x direction is ey, = om = 0.0015 and in the y direction is ey, = —ve,, = —0.00045. Hence the 
length of the strip increases by 0.0015 x 4 = 0.006 and the width reduces by 0.00045 x 2 = 0.0009. 
The total forces acting on the left hand and right hand edges are both 6 (stress times edge length, the 
strip thickness being treated as unity). 


y A E B 
Poe 


Figure 6.1 Thin strip in simple tension 
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Although the strip is in force equilibrium, the problem cannot yet be solved using boundary elements 
because the boundary conditions are given purely in terms of stresses, with no displacements prescribed. 
This difficulty can be overcome by applying suitable point constraints. The constraints adopted here are 
to fix the midpoint E of side AB in both the x and y directions, while constraining midpoint F of side 
DC not to move in the x direction (but to remain free to move in the y direction). The point constraints 
prevent rigid body movement of the strip as a whole, without giving rise to any point forces at either 


point E or point F 


To obtain a quadratic boundary element solution to the problem, the boundary has first to be divided 
into segments. Four segments along the four edges of the domain are appropriate, because they are 
straight lines between corners and the boundary condition over each one is uniform. The ends of the 
four segments are the corner points labelled A, B, C and D (in that order) in Figure 6.1. The choice of 
starting point is arbitrary, but then the direction of numbering must keep the domain to the left. The 


data file DATA to solve the problem using one element (three nodes) per side of the domain is as follows 


SIMPLE TENSION 


STRESS (Plane stress selected) 
2000. 0.3 (Young’s modulus and Poisson’ ratio) 
1 (Number of boundaries) 

(Number of segments on the boundary) 
Ae Ov Ae 2. Oe 2. (Co-ordinates of the segment ends) 
(Radii of curvature, 
number of elements 


and element length ratio 


PP PRP oO 
7 
PoP PP 


for each of the 4 segments) 
2 (Two prescribed stresses, three point constraints) 
Normal stress = 3 on segment 2) 


Normal stress = 3 on segment 4) 


Node 2 constrained in y direction) 


a NW DNB FSF NY CO GO OO Oo oO Oo LH 


PF NO FP WwW W SN 
fo) 


( 
( 
(Node 2 constrained in x direction) 
( 
(Node 6 constrained in x direction) 


The annotations on the right have been added for explanation, and are not part of the data file. Note that 


the prescribed zero values of stresses (on DC and AB) do not have to be entered explicitly. 
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idaryElementPrograr 


GEOMETRIC DATA FOR THE MESH 


NUMBER OF ELEMENTS = 4 


NUMBER OF NODAL POINTS = 8 


COORDINATES OF THE NODAL POINTS 
ne x Y I x Y 
1 0.0000E+00 0.0000E+00 2 0.2000E+01 0.0000E+00 
3. 0.4000E+01 0.0000E+00 4 0.4000E+01 0.1000E+01 
5 0.4000E+01 0.2000E+01 6 0.2000E+01 0.2000E+01 
7 0.0000E+00 0.2000E+01 8 0.0000E+00 0.1000E+01 
ELEMENT NODE NUMBERS 
M ND1 ND2 ND3 M ND1 ND2 ND3 
1 t 2 3 2 3. 4 5 
3 «5 6 7 4 7 . 1 


With one element per segment, 4 elements and 8 nodes are created. Both elements and nodes are 
numbered as for the segments, anticlockwise from the origin, as shown in Figure 6.2. Element numbers 


are shown circled. Element end nodes are shown as small solid circles, mid-side nodes as open circles. 


Figure 6.2 Discretisation of simple tension test problem 


The RESULTS data file output by the program is 


QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 
SIMPLE TENSION 

UNDER PLANE STRESS CONDITIONS 

YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.300 

PRESCRIBED STRESS BOUNDARY CONDITIONS 


SIGNN = 0.3000E+01 SIGSN = 0.0000E+00 ON ELEMENTS 2 TO 2 


SIGNN = 0.3000E+01 SIGSN = 0.0000E+00 ON ELEMENTS 4 TO 4 


NODE 2 CONSTRAINED WITH U = 0 
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NODE 2 CONSTRAINED WITH V 


0 
NODE 6 CONSTRAINED WITH U = 0 


NODAL POINT DISPLACEMENTS AND TRACTIONS 


I U V TX PY. 

1 -0.299965E-02 -0.101365E-05 -0.300000E+01 0.000000E+00 
2 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 
3 0.299965E-02 -0.101577E-05 0.300000E+01 -0.799680E-06 
4 0.300182E-02 -0.450233E-03 0.300000E+01 0.000000E+00 
5 0.299965E-02 -0.899451E-03 0.150000E+01 -0.399840E-06 
6 0.000000E+00 -0.900459E-03 0.000000E+00 0.000000E+00 
7 -0.299965E-02 -0.899450E-03 -0.300000E+01 0.000000E+00 
8 -0.300182E-02 -0.450231E-03 -0.300000E+01 0.000000E+00 


STRESSES AT THE NODES OF THE ELEMENTS 


M IN SIGNN SIGSS SIGSN SIGE 

1 1 0.000000E+00 0.299965E+01 0.000000E+00 0.299965E+01 
1 2 0.000000E+00 0.299965E+01 0.000000E+00 0.299965E+01 
1 3 0.000000E+00 0.299965E+01 0.000000E+00 0.299965E+01 
2 1 0.300000E+01 0.156770E-02 0.000000E+00 0.299922E+01 
2 2 0.300000E+01 0.156531E-02 0.000000E+00 0.299922E+01 
2 3 0.300000E+01 0.156287E-02 0.000000E+00 0.299922E+01 
3 1 0.000000E+00 0.299965E+01 0.000000E+00 0.299965E+01 
3 2 0.000000E+00 0.299965E+01 0.000000E+00 0.299965E+01 
3 3 0.000000E+00 0.299965E+01 0.000000E+00 0.299965E+01 
4 1 0.300000E+01 0.156229EF-02 0.000000E+00 0.299922E+01 
4 2 0.300000E+01 0.156415E-02 0.000000E+00 0.299922E+01 
4 3 0.300000E+01 0.156595E-02 0.000000E+00 0.299922E+01 


FORCES ACTING ON THE BOUNDARY SEGMENTS 


SEGMENT FX FY 

1 0.00000E+00 0.00000E+00 
2 0.60000E+01 -0.53312E-06 
3 0.00000E+00 0.00000E+00 
4 -0.60000E+01 0.00000E+00 
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The prescribed stress conditions are a normal value of 3 and shear value of zero over element 2 (nodes 3, 
4 and 5) on the right end edge (segment 2), and a normal value of 3 and shear value of zero over element 
4 (nodes 7, 8 and 1) on the left hand edge (segment 4). Segments (and elements) 1 and 3 defining the 
top and bottom edges had no boundary condition prescribed by the DATA file, so zero stress conditions 


are assumed. 


The computed displacements show that at nodes 3, 4 and 5 (on the right hand edge) the displacements 
in the direction are 0.0029996, 0.0030018 and 0.0029996, respectively, while at nodes 7, 8 and 1 (on the 
left hand edge) they are -0.0029996, -0.0030018 and -0.0029996. The extensions of the strip between 
node pairs 3 and 1, 4 and 8 and 5 and 7 are therefore 0.0059992, 0.0060036 and 0.0059992, respectively. 
These compare very favourably with the exact value of 0.006. Note the slight differences between element 
end nodes and element midside nodes, which were also experienced in the quadratic element analysis of 
potential problems. The displacement in the y direction of node 6 (which is vertically above constrained 
node 2) is -0.00090046, which compares with the exact value of -0.0009. 


The value of Ox, (SIGSS) computed at all the nodal points on the top and bottom edges is 2.9996, 
compared with the exact value of 3. Also, small tensile values of oy, (SIGSS) of about 0.00156 are 
computed at all the nodes on the right and left hand side edges, where they should be zero. The forces on 
the left and right hand edges are exact at 6.0000, which is to be expected because they were prescribed 


in the boundary conditions. 
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The results are highly accurate for this very simple problem. Using just one quadratic element and three 


nodes on each of the four edges of the boundary gives results which are almost correct to four significant 


figures. For this very simple problem it is not necessary to refine the mesh because the accuracy of the 


results is already adequate for most practical engineering purposes. 


Using displacement boundary conditions 


An alternative approach to the simple tension problem shown in Figure 6.1 is to apply displacement 


boundary conditions to one end of the strip, say edge AD. Unfortunately, this normally changes the 


nature of the problem because contraction along AD is prevented. The only exception is if the value of 


Poisson’ ratio is set to zero, rendering the problem one-dimensional. 


The RESULTS data file output for this case is 


PLE 


ENSION 


ER PLANE 


ST 


RE 


DULUS 


DP 


ED DISE 


DRATIC BOUNDARY ELE 


CONDITIONS 
0.2000E+04 


EM 


OE+00 


ED STR 


NO 


Oo WY WD oO FPF WB NY FP WH 


NM ND F&F 


U 


0.000000E+00 
.300082E-02 


[> ame = aa <3 > a > a <P 


-600076E 
-600260E 
-600076E 
-300083E 


-000000E+00 
-000000E+00 


ES AT THE 


NO FR WwW 


ESS 


-3000E+01 SI 


0.0000E+00 ON 


GSN 


EM 


Vv 


-000000E+00 
.376321E-07 
-847336E-06 
-344766E-08 
-840583E-06 
-407261 
-000000E+00 
-000000E+00 


E-0O7 


-0.300536 
0.000000] 
0.300000] 


-0.300536 
-0.299728 


POISSON'S RATIO 


ELEM 


E+00 


BOUNDARY CONDITIONS 
0.0000! 
ENTS AND TRACTIONS 


TX 


= 
G 
py 


ON 


ENT SOLUTION FOR TWO DIM 


ENT BOUNDARY CONDITIONS 


ENSIONAL ELASTIC PROBLE 


ENTS 4 TO 4 


ENTS 2 TO 2 


E+01 


+00 


E+01 


0.300000E+01 


0.150000! 


E+01 


0.000000E+00 


NODES OF THE 


SIGNN 


000000E+00 
-000000E+00 
-000000E+00 
-300000E+01 
-300000E+01 


SIGSS 
-300127E+01 
.300038E+01 
-299949E+01 
.168764E-02 
.168792E-02 
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ro 


oOo Oo Oo Oo © 


0000001 
. 0000001 
0000001 
0000001 
. 0000001 


ENTS 


SIGSN 


E+01 
E+01 


E+00 
E+00 
E+00 
E+00 


E+00 


TY 
0.254296E-03 
0.000000E+00 

-0.799680E-06 
0.000000E+00 
-0.399840E-06 
0.000000E+00 
-0.251318E-03 
-0.181283E-06 


SIGE 
0.300127E+01 
0.300038E+01 
0.299949E+01 
0.299916E+01 
0.299916E+01 


0.000 


2 3 0.300000E+01 0.168820E-02 0.Q000000E+00 0.299916E+01 
3 1 0.000000E+00 0.299949E+01 0.000000E+00 0.299949E+01 
3 2 0.000000E+00 0.300038E+01 0.000000E+00 0.300038E+01 
3 3 0.000000E+00 0.300127E+01 0.000000E+00 0.300127E+01 
4 1 0.300536E+01 0.000000E+00 0.251318E-03 0.300536E+01 
4 2 0.299728E+01 0.000000E+00 0.181283E-06 0.299728E+01 
4 3 0.300536E+01 0.000000E+00 -0.254296E-03 0.300536E+01 


FORCES ACTING ON THI 


BOUNDARY SEGMENTS 


Gl 


SEGMENT FX FY 
1 0.00000E+00 0.00000E+00 
2 0.60000E+01 -0.53312E-06 
2 0.00000E+00 0.00000E+00 
4 -0.60000E+01 0.75079E-06 


The extensions of the strip computed at nodes 3, 4 and 5 on edge BC are 0.0060007, 0.0060026 and 
0.0060007 (exact 0.006), the o,, stresses at nodes on the top and bottom edges of the strip are 3.0013, 
3.0004 and 2.9995 (exact 3), and the force on edge AD is exact at 6.0000. The accuracy remains good 


when displacements are prescribed rather than stresses. 


6.2.2 Pure shear 


In order to explore shear loading, Figure 6.3 shows the cross-section of a rectangular block of material 
which is long in the direction normal to the cross-section. The material concerned again has a Young's 
modulus of 2000 and a Poisson's ratio of 0.3. The base of the block, AB, is rigidly constrained with 
prescribed zero values of displacements. The other edges have shear stresses of magnitude 3 applied in 
such a way as to give rise to a state of pure shear. In other words, on sides BC and DA o,, = +3, and 
on side CD o,,, = —3, bearing in mind that co-ordinate S is measured along the boundary keeping the 


domain always to the left. 


Figure 6.3 Block under pure shear 
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The shear modulus for the material is 


E 2000 


@=— Se 692 
daw) Ces 


and the shear strain due to a shear stress of 3 units is 3/769.23 = 0.0039. With a block height of 2, the 
displacement of top surface DC is u = 2 x 0.0039 = 0.0078. 


The RESULTS data file output for this case is 


QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 


PURE SHEAR 


UNDER PLANE STRAIN CONDITIONS 


YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.300 


PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS 
U = 0.0000E+00 V = 0.0000E+00 ON ELEMENTS 1 TO 1 


PRESCRIBED STRESS BOUNDARY CONDITIONS 


SIGNN = 0.0000E+00 SIGSN = 0.3000E+01 ON ELEMENTS 2 TO 2 
SIGNN = 0.0000E+00 SIGSN = -0.3000E+01 ON ELEMENTS 3 TO 3 
SIGNN = 0.0000E+00 SIGSN = 0.3000E+01 ON ELEMENTS 4 TO 4 


> Apply now 


REDEFINE YOUR FUTURE 
AXA GLOBAL GRADUATE 
PROGRAM 2015 


© Photononstop 


redefining / standards 


Download free eBooks at bookboon.com 


108 Click on the ad to read more 
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NODAL POINT DISPLACEMENTS AND TRACTIONS 
I U V TX TY 
1 0.000000E+00 0.000000E+00 -0.300535E+01 -0.184857E-02 
2 0.000000E+00 Q.000000E+00 -0.299735E+01 0.206969E-06 
3 0.000000E+00 0.000000E+00 -0.300535E+01 0.184803E-02 
4 0.390155E-02 0.144327E-05 0.000000E+00 0.300000E+01 
5 0.780091E-02 0.117580E-06 0.150000E+01 0.150000E+01 
6 0.780417E-02 0.283357E-09 0.300000E+01 0.000000E+00 
7 0.780091E-02 -0.118911E-06 0.150000E+01 -0.150000E+01 
8 0.390155E-02 -0.144314E-05 0.000000E+00 -0.300000E+01 
STRESSES AT THE NODES OF THE ELEMENTS 
M IN SIGNN SIGSS SIGSN SIGE 
1 1 0.184857E-02 0.792242E-03 -0.300535E+01 0.520542E+01 
1 2 -0.206969E-06 -0.887008E-07 -0.299735E+01 0.519156E+01 
1 3 -0.184803E-02 -0.792012E-03 -0.300535E+01 0.520542E+01 
2 1 0.000000E+00 0.621710E-02 0.300000E+01 0Q.519616E+01 
2 2 0.000000E+00 0.129209E-03 0.300000E+01 0.519615E+01 
2 3  0.000000E+00 -0.595411E-02 0.300000E+01 0.519616E+01 
3 1 0.000000E+00 -0.716504E-02 -0.300000E+01 0.519616E+01 
3 2 0.000000E+00 0.000000E+00 -0.300000E+01 0.519615E+01 
3.3 0.000000E+00 0.716504E-02 -0.300000E+01 0.519616E+01 
4 1 0.000000E+00 0.595146E-02 0.300000E+01 0.519616E+01 
4 2 0.000000E+00 -0.130672E-03 0.300000E+01 0.519615E+01 
4 3  0.000000E+00 -0.621280E-02 0.300000E+01 0.519616E+01 


FORCES ACTING ON THE BOUNDARY SEGMENTS 


SEGMENT FX FY 
1 -0.12000E+02 0.19397E-06 
2 0.53312E-06 0.60000E+01 
3 0.12000E+02 -0.26656E-06 
4 0.Q00000E+00 -0.60000E+01 


The computed displacements at nodes 5, 6 and 7 on the top surface DC are 0.0078009, 0.0078042 and 
0.0078009 (exact 0.0078). The computed shear stresses at nodes 1, 2 and 3 on the bottom surface AB 
are -3.00535, -2.9973 and -3.0053 (exact 3), and the total shear force on this surface is exact at -12.000. 


Download free eBooks at bookboon.com 


109 


It is worth noting that had the shear stresses on sides BC and DA been omitted in Figure 6.3, although the 
block might appear to be simply sheared, the actual state of stress would have been much more complex. 
In particular, at corners A and B there is a conflict between finite shear stress along the constrained base 


and zero shear stress on the side of the block. Equation 1.1 requires shear stresses to be complementary. 


6.2.3 Steel component 


Figure 6.4 shows a thin steel component with elastic properties of Young’s modulus 207 GN/m’, Poisson's 
ratio 0.30. The bottom edge is subject to a uniform normal stress of 100 MN/m/”, and the top edges 
to uniform normal stresses of 75 MN/m? (giving force equilibrium in the vertical direction). The top 
edges are also subject to outwardly directed shear stresses of 30 MN/m’, giving force equilibrium in the 


horizontal direction. 


400 


250 


All dimensions in mm 


Figure 6.4 A steel component under load 


This problem involves somewhat more complex geometry and varying stresses, and is therefore likely to 
require rather more mesh refinement than the very simple problems treated so far. The object is to find 
the direct stress in the surface at point E (mid way between points D and F on the semi-circular arc), 
and the maximum von Mises equivalent stress anywhere on the boundary (and hence anywhere in the 


component). There is no exact analytical solution for comparison. 


The boundary of the domain is conveniently divided into six segments: AB, BC, CD, DEF, FG and GA. 
For convenience, the same number of boundary elements are used on each segment. Displacement 
constraints are required to prevent rigid body movement: for example, point A fixed in both directions, 
and point B constrained not the move in the direction. Table 6.1 shows the effects of varying the number 
of nodes per segment from 4 to 40. Less than 4 gives poor geometric representation of the semi-circle 
DEF (with each element then representing an arc of more than 45°). 
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Elements Os; at E max.d, 


per segment MN/m? MN/m? 
4 -111.3 152.2 
(140.2) 

6 -119.9 151.5 
(148.7) 

10 -123.8 153.7 
(153.7) 

20 -125.1 152.9 
(152.9) 

40 -125.3 152.7 
(152.7) 


Table 6.1 Stress results for varying levels of mesh refinement 


In each case point E, which is on the vertical line of symmetry for the problem, is located at the common 
node of two elements, and the value of %s; obtained from each element is the same (to at least four 
significant figures). On the other hand, the maximum value of equivalent stress, which occurs in the 
regions of both points H and I in Figure 6.4, does not necessarily occur at a node at all. This means that 
a reasonable degree of mesh refinement is needed to capture a good approximation to the maximum. In 
each case the maximum is detected at an element end node, and the two values of equivalent stress from 
the two elements meeting at that node are not necessarily the same. In the table the values in brackets 


are those from the adjacent element. 


There are two criteria to keep in mind when testing for convergence with mesh refinement. Firstly, 
whether a value of stress changes with increasing number of elements. Secondly, whether stress values 
at an element end node (in this case the maximum equivalent stress) are the same for both elements 
meeting at that node. With these in mind, the second criterion in this problem is satisfied at 10 elements 
per segment or more, whereas the first is only satisfied (to three significant figures) at 20 elements per 
segment or more. The computed direct stress in the surface at point E is 125 MN/m? (compressive) and 
the maximum von Mises equivalent stress 153 MN/m* (actually the tensile direct stresses in the surfaces 


in the regions of points H and I). 


6.2.4 A thick-walled cylinder test problem 


The next test is designed to demonstrate the ability of the program to solve problems with curved 
boundaries and those with more than one boundary. Figure 6.5 shows the cross-section of a long thick- 
walled circular cylinder with internal and external radii of % = 4 and 12 = 8. The cylinder is subject to 
an internal pressure of p = 3 and no pressure on the outside. The Young’s modulus and Poisson's ratio 
of the material of the cylinder are 2000 and 0.3, respectively. Cylinders are widely used in engineering 
practice to contain fluids under pressure: from pipes for relatively low pressures to thick-walled vessels 


for high-pressure chemical reactors. 
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Figure 6.5 Thick-walled cylinder problem 


The analytical solution to the thick-walled cylinder problem, under either plane stress or plane strain 


conditions, gives a general result for the radial and hoop stresses as 


B B 
Orr = =. 069 = A+ 2. (6.5) 
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where r is radial distance from the centre of the cylinder, and and are constants. Stress 0; is the direct 
radial stress in the r direction while hoop stress %@ is, as its name suggests, the direct stress in the 
circumferential direction at any given radial position (for example, along the inner or outer surfaces of 
the cylinder). The values of the constants depend on the boundary conditions applied to the cylinder. 


In the present problem, these are 


Orn =—p at r=n7 and 6 =0 at r=% (6.6) 
2 
_ p : —_ Pl a2 
Hence A= aD: and B iK2—-1 where K 7 (6.7) 


and the hoop stresses at the inner and outer cylinder surfaces are 


(K* +1) 2p 
099 = aes at r=" and 099 = (2-1) at r=" (6.8) 


In this problem K = 2 and p = 3, and these hoop stresses are 5 and 2, respectively. 


These results for stresses do not depend on the values of either Young’s modulus or Poisson's ratio. Radial 


displacement, u,, can be found from the hoop strain 


C99 = 7 = = (op —v"o,) (6.9) 

At the inner cylinder surface u, = a (5 + 3v*) (6.10) 

and at the outer surface U, = a (6.11) 
From Equations 5.7 

"= ose = 219780 and v= + = 0.428571 (6.12) 


and the displacements at the inner and outer cylinder surfaces are 


uy, = 0.011440 and uy = 0.007280 
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Bearing in mind that program BEM2EQ cannot accept boundary segment specifications with angles 
greater than , the simplest mesh to represent the geometry shown in Figure 6.5 requires four segments, 
two on the outer boundary with end points arbitrarily chosen at (8, 0) and (-8, 0), and two on the inner 
boundary with end points at (4, 0) and (-4, 0). These end points are the points F, G, H and I in Figure 
6.5. The origin for the global Cartesian co-ordinates is the centre of the cylinder. The following is an 
annotated DATA input file for this problem, for the case of four elements (9 nodes) per semi-circular 
segment, two per quadrant. This means that each element represents an arc of by means of a parabolic 


(quadratic) curve passing through both the ends and midpoint of the arc. 


THICK-WALLED CYLINDER 


STRAIN (Plane strain selected) 

2000. 0.3 (Young’s modulus and Poisson’s ratio) 

2 (Number of boundaries) 

2 (Number of segments on first boundary) 

8. 0. -8. 0. (Co-ordinates of the segment ends) 

Sie, A ses (Radii of curvature, numbers of elements 

8. 41. and length ratios) 

2 (Number of segments on second boundary) 
4. 0. -4. 0. (Co-ordinates of the segment ends) 

-4, 41. (Radii of curvature, numbers of elements 
=4; 4 1. and length ratios) 

0 2 3 (2 prescribed stresses and 3 point constraints) 
3-3. 0. (Normal stress = -3 on segment 3) 

A. = 3:3, 05 (Normal stress = -3 on segment 4) 

11 (Node 1, point F, constrained in x direction) 
12 (Node 1, point F, constrained in y direction) 
9.2 (Node 9, point G, constrained in y direction) 


As in the case of the simple tension problem of Section 6.2.1, sufficient displacements constraints must 
be prescribed to prevent rigid body movement of the whole cylinder. Those chosen are no movement 
in either the or direction for point F, and no movement in the direction for point G. In terms of node 


numbers, for the present mesh of 4 elements per segment these are 1 and 9. 
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Table 6.2 shows the results obtained for meshes ranging from one element per quadrant up to four 
elements per quadrant. For the computed stresses, two values are shown, the first for element end nodes, 
and the second for midside nodes. The computed displacements are for the highest points on surfaces, 


J and K in the figure, in the direction. 


Elements At inner surface At outer surface 
per quadrant ory) Uy lorry Uy 

Exact 5 0.011440 2 0.007280 

1 4.8985 0.011029 1.9653 0.007030 
4.7743 1.9313 

2 4.9940 0.011417 1.9967 0.007213 
4.9875 1.9965 

3 5.0006 0.011437 1.9996 0.007279 
4.9984 1.9996 

4 5.0014 0.011440 2.0001 0.007280 
5.0001 2.0001 


Table 6.2 Results for the thick-walled cylinder problem 
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Since the problem is axi-symmetric it is to be expected that the computed stresses at all the nodes on the 
inner boundary would be identical, and similarly for all the nodes on the outer boundary. This proves to 
be the case, but for expected fluctuations between the values at element end nodes and midside nodes. 
At four elements per quadrant the results are almost correct to four significant figures, which is adequate 
for most practical purposes. In this case each element covers a circular arc of 22.5°, which demonstrates 
the ability of quadratic elements to accurately model curved boundaries. Thirty two elements having 64 


nodes are sufficient to model this multiply-connected solution domain. 


6.3 An Example: Confined Compression of a Rubber Block 


A practical problem for solution by program BEM2EQ is the compression of a rubber block which is 
confined between rigid surfaces. It has important engineering applications in the bearings for bridges. 
Bridges are often supported on rubber pads. These are relatively flexible in shear, allowing the bridge to 
expand or contact in length with changes in ambient temperature and the load it is carrying. But they 
are much stiffer in compression. The reason for this is that rubber, which is much more flexible than 
metals, is incompressible (or very nearly incompressible), having a Poisson's ratio of 0.5. It is worth 
noting that many finite element methods are incapable of analysing incompressible material problems 


under plane strain conditions. 
D 
ny 
x 
A E 


Figure 6.6 Rubber block between rigid surfaces 


C 
B 


Figure 6.6 shows a rubber block sandwiched between rigid surfaces, to which it is bonded. The height 
and width of the block are H and W, respectively. Consider the case where W/H = 10.. The dimension 
of the block in the direction normal to the plane shown is large, so that a state of plane strain may be 
assumed. With the bottom surface AB fixed in both global co-ordinate directions, the top surface is to 
be lowered by 0.1% of the thickness H. The (maximum) compressive stress at point E is to be computed, 
together with the total compressive force on the surface AB, both to be expressed as multiples of the 


stress and force which would exist if the block were in a state of simple compression. 


If the Young’s modulus is taken as 1000, a 0.1% simple compressive strain would give a stress of 
0.001 x 1000 = 1. If block width is taken as 1, then the total forces on the upper and lower surfaces 
would also be 1 under the same conditions. So in the boundary element model it is convenient to 
take W = 1,H = 0.1 and an upper surface vertical displacement of vy = —0.0001 (0.1% of 0.1). The 


computed stresses and forces will be the required multiples. 
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The boundary of the domain is conveniently divided into four segments: AB, BC, CD and DA. With 
a wide disparity in the lengths of the segments, it is no longer appropriate to use the same number of 
elements on each segment. Best results are normally obtained when the lengths of adjacent elements are 
reasonably similar (not differing by more than a factor of about 2 or 3). In this problem it is appropriate 
to have ten times as many elements on AB and CD as on BC and DA. This means that the coarsest mesh 
which can be used has 10 elements each on AB and CD, and only one each on BC and DA. This is also 
desirable from the point of view that the domain is slender and it is important that opposite sides are 
not too close together in relation to element sizes (for the reasons discussed in Section 4.2.1 for potential 
problems). With 10 elements on each of the horizontal edges, their distance apart is equal to the common 
lengths of the elements. Mesh refinement (keeping the 10:1 element number ratio between the horizontal 


and vertical segments) will only improve the situation. 


Table 6.3 shows the effects of varying the number of elements on the horizontal segments from 10 to 50 
(with the numbers of elements on the vertical segments always one tenth of these). 


Elements Oyy atE Total force 
per segment on AB 
10 -50.43 -34.51 
20 -50.15 -34.03 
30 -50.18 -33.95 
40 -50.19 -33.92 
50 -50.19 -33.91 


Table 6.3 Stress and force results for varying levels of mesh refinement 


There is no change in the maximum stress of 50.2 beyond 20 elements per horizontal segment and 
the total force has settled to 33.9 by 30 elements per segment. As multiples of the values expected for 
simple compression these are large numbers, and would have been much higher if larger ratios had 
been considered. The explanation is that, because the rubber material is assumed to be completely 
incompressible, the reduction in the height of the block caused by bringing the confining surfaces closer 
together can only be accommodated by material being forced out horizontally at the side edges, BC 
and DA. Because there are no displacements at the corners A, B, C and D the rubber is forced out into 
curved bulges. In this case the maximum horizontal displacements at the centres of BC and DA are each 
computed as 0.00072, 7.2 times the relative displacement of the confining surfaces. The driving forces 
required to produce this deformation are provided by the stress gradients in the horizontal direction 
moving away from the centre of the block. From a maximum of 50.2 at point E, both o,, and o,, (the 
state of stress in the block is hydrostatic, with these stresses equal) follow a roughly parabolic horizontal 
distribution to finish at close to one at sides BC and DA. The average stress is 33.9, numerically equal 


to the total applied force. 
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6.4 An Example: Stress Concentration at a Hole in a Flat Plate 


As a final example, consider the classic problem of the stress concentration created at a small circular 
hole at the centre of a thin flat square plate under tension. This is often used as a test of any numerical 
method of stress analysis. Figure 6.7 shows the geometry and loading. The radius of the hole is , and the 
width and height of the plate are 2: the ratio of hole diameter to plate width is . For present purposes this 
ratio is taken as 1/20. A uniform tensile stress is applied to both the top and bottom edges of the plate. 


r a 4 


o 


Figure 6.7 Square plate with a small central hole 


There is no exact analytical solution for the case of a plate of finite width and height, although one does 
exist for an infinitely large plate. This gives a maximum (tensile) hoop stress at both points E and G of 
30, together with a (compressive) hoop stress at both points F and H of -o. It also indicates that there 
are rapid changes in stresses and deformations close to the hole, which make this a relatively challenging 
application for any numerical method. The present problem is for a plate of finite width, although the 
width to diameter ratio is relatively high. The maximum stress concentration factor (ratio of local stress 


to remote applied stress) can be expected to be a little higher than 3. 


The outer boundary of the domain is conveniently divided into four segments: AB, BC, CD and DA, and 
the inner boundary into two: EHG and GFE. On the outer boundary the stresses and displacements do 
not vary much and only a few boundary elements are required: four per segment is sufficient. Rather 
more elements are likely to be required on the inner boundary. Point constraints are required to prevent 
rigid body movement: point A is constrained in both directions, while point B is constrained not to 


move in the y direction. 
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The DATA file for this problem with four elements per semi-circular segment (two per quadrant) is 


STRESS CONCENTRATION 
STRESS 
2000. 0.3 


=0'9 025 05 O25 O25 0.5. 0.25 
Les 


oO Go & oS 
PP PF BP OF 


1. 
1. 
1 
2 
0.025 0. -0.025 0. 
=0.025 4 1. 
=0.,025 4, 1. 

3 


QO. 
QO. 


Oo FP FPF WF Oo 
NO NM FP FP FP ND 


(Plane stress selected) 

(Young’s modulus and Poisson's ratio) 
(Number of boundaries) 

(Number of segments on first boundary) 
(Co-ordinates of the segment ends) 
(Radii of curvature, numbers of elements 


and length ratios) 


(Number of segments on second boundary) 
(Co-ordinates of the segment ends) 

(Radii of curvature, numbers of elements 

and length ratios) 

(2 prescribed stresses and 3 point constraints) 
(Normal stress = 1 on segment 1) 

(Normal stress = 1 on segment 3) 

(Node 1 constrained in x direction) 


(Node 1 constrained in y direction) 


(Node 9 constrained in y direction) 
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The overall dimensions of the plate are taken as 1, so that a=0.025 and b=0.5. The Young’s modulus and 
Poisson's ratio are arbitrarily chosen as 2000 and 0.3, respectively, although they do not affect the results 


for stress concentration. 


An edited version of the RESULTS data file output produced is 


QUADRATIC BOUNDARY ELEMENT SOLUTION FOR TWO DIMENSIONAL ELASTIC PROBLEM 
STRESS CONCENTRATION 

UNDER PLANE STRESS CONDITIONS 

YOUNG'S MODULUS = 0.2000E+04 POISSON'S RATIO = 0.300 

PRESCRIBED STRESS BOUNDARY CONDITIONS 


SIGNN 


0.1000E+01 SIGSN 


0.0000E+00 ON ELEMENTS 1 TO 4 


SIGNN 


0.1000E+01 SIGSN 


0.0000E+00 ON ELEMENTS 9 TO 12 
0 
D 0 
NODE 9 CONSTRAINED WITH V = 0 
R 


NODE 1 CONSTRAINED WITH U 


NODE 1 CONSTRAINED WITH V 


STRESSES AT THE NODES OF THE ELEMENTS 


M IN SIGNN SIGSS SIGSN SIGE 

17 1 0.000000E+00 0.301684E+01 0.000000E+00 0.301684E+01 
17 2 0.000000E+00 0.242750E+01 0.000000E+00 0.242750E+01 
17 3 0.000000E+00 0.105800E+01 0.000000E+00 0.105800E+01 
18 1 0.000000E+00 0.945216E+00 0.000000E+00 0.945216E+00 
18 2 0.000000E+00 -0.424362E+00 0.000000E+00 0.424362E+00 
18 3 0.000000E+00 -0.101370E+01 0.000000E+00 0.101370E+01 
19 1 0.000000E+00 -0.101370E+01 0.000000E+00 0.101370E+01 
19 2 0.000000E+00 -0.424375E+00 0.000000E+00 0.424375E+00 
19 3 0.000000E+00 0.945128E+00 0.000000E+00 0.945128E+00 
20 1 0.000000E+00 0.105804E+01 0.000000E+00 0.105804E+01 
20 2 0.000000E+00 0.242751E+01 0.000000E+00 0.242751E+01 
20 3 0.000000E+00 0.301677E+01 0.000000E+00 0.301677E+01 
21 1 0.000000E+00 0.301691E+01 0.000000E+00 0.301691E+01 
21 2 0.000000E+00 0.242751E+01 0.000000E+00 0.242751E+01 
21 3 0.000000E+00 0.105795E+01 0.000000E+00 0.105795E+01 
22 1 0.000000E+00 0.945229E+00 0.000000E+00 0.945229E+00 
22 2 0.000000E+00 -0.424373E+00 0.000000E+00 0.424373E+00 
22 3 0.000000E+00 -0.101371E+01 0.000000E+00 0.101371E+01 
23 1 0.000000E+00 -0.101371E+01 0.000000E+00 0.101371E+01 
23 2 0.000000E+00 -0.424386E+00 0.000000E+00 0.424386E+00 
23 3 0.000000E+00 0.945244E+00 0.000000E+00 0.945244E+00 
24 1 0.000000E+00 0.105804E+01 0.000000E+00 0.105804E+01 
24 2 0.000000E+00 0.242751E+01 0.000000E+00 0.242751E+01 
24 3 0.000000E+00 0.301683E+01 0.000000E+00 0.301683E+01 


FORCES ACTING ON THE BOUNDARY SEGMENTS 
Download free eBooks at bookboon.com 


120 


Boundary Element Methods for Engineers: 


Part II: Plane Elastic Problems QuadraticBoundaryElement Program forPlaneElasticProblems 
SEGMENT FX BY 

1 -0.56196E-07 -0.10000E+01 

2 0.00000E+00 0.00000E+00 

3 0.56196E-07 0.10000E+01 

4 0.00000E+00 0.00000E+00 

5 0.00000E+00 0.00000E+00 

6 0.00000E+00 0.00000E+00 


Results for nodal point displacements and tractions have been omitted, together with stresses on the 
outer boundary (elements 1 to 16). The hoop stress (SIGSS) at point E is either that for the first node of 
element 17 (3.01684) or the third node of element 24 (3.01683). At point G it is either that for the third 
node of element 20 (3.01677) or the first node of element 21 (3.01691). The four values are consistent, 
and give an average to five significant figures of 3.0168. The hoop stress at point F is either that for the 
third node of element 18 (-1.01370) or the first node of element 19 (-1.01370). At point H it is either 
that for the third node of element 22 (-1.01371) or the first node of element 23 (-1.01371). The four 
values are consistent, and give an average to five significant figures of -1.01371. With only two elements 


per quadrant of the hole in the plate the stress values are very close to what is expected. 
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Elements Oyy Oe 
per quadrant at E and G at F and H 
2 3.0168 -1.0137 
3 3.0206 -1.0163 
4 3.0211 -1.0167 
6 3.0215 -1.0168 


Table 6.4 Stress concentration results for varying levels of mesh refinement at the hole 


Table 6.4 shows the effects of refining the mesh around the hole. Both stress values increase in magnitude 
a little and converge to 3.021 at the sides of the hole and -1.017 at the top and bottom of the hole. 
Converged results are achieved with only a modest number of boundary elements: at 6 elements per 
quadrant the total used is 40. With rapid variations of stresses and deformations near the hole, this is the 
type of problem where the ability to vary the element distributions along segments (to give small elements 
at the stress concentrations) might have been expected to offer significant benefits. The quality of the 
results from using elements of constant size suggests, however, that this is not necessary. The situation 
would have been different if symmetry had been invoked to model only one quarter of the plate (say, the 
first quadrant in the plane). Element size variation on the lines of symmetry would have been required 


to accommodate both small elements near the hole and much larger elements at the outer boundary. 


6.5 Discussion 


It is clear that quadratic boundary elements applied to elastic stress analysis problems offer good 
accuracy with relatively small numbers of elements. The scope of the program described in this chapter 
is deliberately limited, to keep the coding reasonably straightforward to understand. Many other features 
can be incorporated: for example, additional types of boundary conditions, displacement and stress 


calculations at internal points, and inclusion of body forces such as gravity and centrifugal effects. 


Another useful facility is the ability is to be able to divide a domain into several distinct regions, with 
elements along the boundaries between them. The different regions can be of different materials, having 
different elastic properties. This makes possible the treatment of not only multi-material problems, but 
also multiple bodies in contact with each other. Another class of problems that can be catered for is that 
of crack problems in linear elastic fracture mechanics. One of the boundary element end nodes can be 
located at the crack tip and the midside nodes in the adjacent elements shifted to force the required stress 
singularity at the crack tip (just as is done in finite element methods). If both faces of the crack have to 
be modelled, then the multi-region facility is very useful. This is because coincident nodal points in a 
single boundary element region are not permitted, and nearly coincident nodes give rise to inaccurate 
analysis. Coincident and nearly coincident nodes in different regions, on the two sides of the crack, 


present no difficulties. 
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Problems 


6.1 


6.2 


6.3 


In the stress concentration problem illustrated in Figure 6.7 (with a/b = 1/20), instead of the 
tensile loading applied to the top and bottom edges of the plate, the hole itself is subjected to 
a uniform internal pressure. Use program BEM2EQ to find the maximum stress concentration 


factor (relative to the applied pressure). Compare this value with the available analytical solutions. 


Consider the long thick-walled cylinder whose cross-section is shown in Figure 6.5, but with its 
outer surface rigidly contained. An internal pressure of p = 1 is applied. The Young’s modulus 
and Poisson's ratio of the material are 1000 and 0.3, respectively. Use program BEM2EQ to find 
both the stresses in the cylinder at the inner and outer surfaces, and the radial displacement of 


the inner surface, and compare with analytical results. 


Figure 6.8 shows a thin rectangular plate with a small semi-circular notch at the centre of one 
horizontal edge. The radius of the notch is 1/20" of the vertical height of the plate. The plate 
is subjected to uniform tensile stresses on its vertical edges. Use program BEM2EQ to find the 


maximum stress concentration at the notch. Take Poisson's ratio to be 0.4. 
r 


0.05R 


— 


Figure 6.8 Notched plate in tension 
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6.4 Figure 6.9 shows a thin square plate with a square hole at its centre. The corners of the hole are 
rounded, with radii of curvature 10% of the hole dimensions, 5% of the overall plate dimensions. 
Use program BEM2EQ to find the maximum stress concentration at the hole. Take Poisson’s 
ratio to be 0.25. 


Figure 6.9 Square plate with a square central hole 


6.5 A cast iron pedestal is trapezoidal in cross-section, with a horizontal base 0.3 m wide, horizontal 
top 0.1 m wide and height 0.1 m, and is symmetrical about the vertical centre line. Plane strain 
conditions may be assumed. The Young’s modulus is 70 GN/m’ and Poisson's ratio 0.25. The top 
surface of the pedestal is uniformly loaded in compression and the base rests on a rigid surface. 
The purpose of the pedestal is to spread the load over the supporting surface. Use program 
BEM2EQ to find the ratio between the maximum compressive stress applied to the supporting 
surface and the stress applied to the top of the pedestal. The surface is rough enough to prevent 


any slipping. Is a coefficient of friction of 0.5 large enough to achieve this? 
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6.6 Figure 6.10 shows a steel tension test piece, which is 3 mm thick. The central section is 20 
mm wide, increasing to 30 mm at the ends where it is gripped. The transition in width is via 
quadrant fillets with radii of 5 mm. Young’s modulus and Poisson's ratio for steel may be taken 
as 207 GN/m? and 0.3, respectively. The ends of the test piece are subjected to uniform stresses 
of 60 MN/m’. Use program BEM2EQ to find the stiffness of the test piece and the maximum 
stress induced. Find also the maximum and minimum axial stresses over the central 30 mm, 
the region where an extensometer is to be attached and where it is intended that the axial stress 


is uniform. 


20 20 
a Le 
~ 100 > 


All dimensions in mm 


Figure 6.10 Tension test piece 
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6.7 


6.8 


6.9 


6.10 


If the edges of the plate in Figure 6.7 (with a/b = 1/20) are subject to pure shear rather than 
simple tension, use program BEM2EQ to find the maximum stress concentration at the hole, 


expressed as a multiple of the applied shear stress. 


Consider the stress concentration problem shown in Figure 6.7 not with a circular hole at the 
centre but with an elliptical hole. The major axis of the ellipse lies along the x axis and is 1/20" 
of the plate width, while the minor lies along the axis and is half the major axis. Use program 


BEM2EQ to find the maximum stress concentration factor. Take Poisson’s ratio to be 0.3. 


Consider the stress concentration problem shown in Figure 6.7 not with a single circular hole at 
the centre but with two holes of equal size (a/b = 1/20), side by side. Their centres are on the x 
axis and equidistant from the centre of the plate. The distance between the centres is four times 
their common radius. Use program BEM2EQ to find the stress concentration factors at the sides 
of the holes. Take Poisson’s ratio to be 0.3. Consider both (a) when the applied uniform tensile 
stresses are in the y direction as shown, at right angles to the line of hole centres, and (b) when 


they are in the x direction, parallel to the line of centres. 


Figure 6.11 shows the cross-section of a long circular steel cylinder, internal and external radii 
100 and 250 mm, respectively. Running the full length of the cylinder is a straight groove which 
is semi-circular in shape with a radius of 5 mm. Young’s modulus and Poisson’s ratio for steel 
may be taken as 207 GN/m* and 0.3, respectively. Use program BEM2EQ to find the maximum 


von Mises equivalent stress in the cylinder when an internal pressure of 120 MN/m’ is applied. 


All dimensions in mm 


TN 
Kae) 


groove 
5R 


Figure 6.11 Thick-walled cylinder with a groove 
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6.11 ‘The long thick-walled internally pressurised cylinder described in Section 6.2.4 has been 
incorrectly manufactured so that it is slightly elliptical in cross-section, with major axis 1% 
greater than minor axis. Use program BEM2EQ to find the ratio between the maximum hoop 
stress (at the internal surface) and the applied pressure, and compare this to the value for the 


intended circular cylinder. 


6.12 Part of the suspension for a motor vehicle involves a torsional spring formed by a circular 
cylindrical rubber bush which is 75 mm long and has inner and outer diameters of 20 mm and 
50 mm, respectively. The outer surface is bonded to a rigid housing, and the inner surface to a 
metal shaft. The rubber has a Young’s modulus of 80 MN/m? and a Poisson's ratio of 0.5. Use 
program BEM2EQ to find the torque which must be applied to the shaft to rotate it by 0.1 radians. 
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7 Further Applications 


In Chapter 1 many continuum mechanics problems in a wide range of engineering subjects were 
categorised according to their governing differential equations into either the harmonic or biharmonic 
types. Problems governed by harmonic equations are, in engineering terms, more often referred to as 
potential problems, such as heat conduction, ideal fluid flow, porous medium flow, torsion and at least 
one form of fluid flow. Part I of the book was devoted to such potential problems, and boundary element 
methods were developed. Part II has been concerned with biharmonic problems in elastic stress analysis, 
particularly two-dimensional plane strain and plane stress. The purpose of this final chapter is to review 


briefly some further applications of boundary element methods. 


7.1 Axi-symmetric Problems 


There are many problems in engineering, both harmonic and biharmonic, which are symmetrical in terms 
of both geometry and boundary conditions about some axis. For example, the thick-walled cylinder test 
problems considered in Sections 3.2.2, 4.2.2 and 6.2.4 were symmetrical about the axis of the cylinder. In 
each case a cross-section of the cylinder at right angles to the axis was treated as the solution domain. It 
would also have been possible to treat as the solution domain a radial cross-sectional plane through the 
wall of the cylinder, which also passed through the cylinder axis. Boundary elements could be used to 
model the boundary of the domain. Interestingly, if a solid rather than hollow axi-symmetric component 
is modelled, so that the axis of symmetry forms part of the boundary of the solution domain, then it is 
not necessary to have boundary elements on the axis. The level of mathematics required to formulate 


axi-symmetric analyses is, however, rather more challenging than for two-dimensional problems. 


Ta Higher-Order Boundary Elements 


It is of course possible to use boundary elements with shape functions of orders higher than quadratic. 
For example, some work has been done with cubic elements, including some which provide a smooth 
transition (continuity of slope) between one element and the next. The high quality of the results that 


can be obtained using quadratic elements, however, means that they are generally the elements of choice. 


7.3 Three-dimensional Problems 


Three-dimensional problems are amenable to solution by a boundary element approach. Again, only 
the boundary of the relevant solution domain is modelled. The surface is covered with two-dimensional 
boundary elements. The equivalent in three-dimensional analysis of the constant element for two- 
dimensional problems is the constant triangle or quadrilateral. Over either such element the potential 
and potential gradient normal to the surface (in a potential problem) or the displacements and tractions 
(in a stress analysis problem) are assumed to be constant over the element, and equal to the value at 
the node at the middle of the element. For potential problems there remains only one unknown at each 
node. For stress analysis there are now three, typically either displacements or tractions in each of the 
three global co-ordinate directions, and three equations are required at each node. 
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Progressing from constant to quadratic elements, the latter can again be either triangular or quadrilateral 
in shape, though not necessarily flat or with straight edges (in the same way that quadratic line elements 
are not necessarily straight lines). Isoparametric elements are generally preferred, with the geometric 
shape being defined in terms of quadratic functions of the three spatial co-ordinates. Triangular elements 
typically have six nodes, quadrilateral eight, one at each corner and one at the centre of each side. 
Quadrilateral elements are often used to model most of a surface, with triangles being useful to deal 
with local geometric features, particularly where a mesh refinement is required to suit a particular 
problem. Indeed, there are situations where a mesh refinement cannot be effected without using at least 


a few triangles. 


While three-dimensional meshes of (curved) triangles and quadrilaterals are significantly more difficult 
to create than two-dimensional meshes of line elements, they still only involve the discretisation of the 
boundary rather than the interior of a domain as well. Using boundary elements, the dimensionality 
of the computational problem is reduced by one. Nevertheless, the amount of computation involved in 
a three-dimensional boundary element analysis is very substantial, particularly because the coefficient 


matrix in the final set of algebraic equations is essentially fully populated with non-zero values. 
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74 Non-Linear Problems 


All the problems so far considered are linear in the sense that they involve the solution of sets of linear 
algebraic equations. Non-linearities can be introduced by either geometric, material property or boundary 
condition effects. Geometric non-linearities arise, for example, in problems involving solid media in 
which the strains are sufficiently large to significantly affect the shape of the solution domain. Such 
problems can be solved by an incremental approach in which the non-linear analysis is replaced by a 
series of linear analyses for progressively increasing external loads, after each of which the boundary 
element mesh geometry is recomputed. Non-linear behaviour due to boundary conditions arises in, for 
example, heat transfer problems with radiation boundary conditions, where the rate of heat transfer to 
or from a boundary depends on the fourth power of the absolute temperature there, rather than the first 
power in a convection boundary condition. In elastic stress analysis problems, contact between two or 
more bodies gives rise to non-linear behaviour: increased contact force between them changes the region 


of contact and hence modifies the boundary conditions. Again incremental approaches are possible. 


Examples of material non-linearities include non-Newtonian fluids and non-linearly elastic or elastoplastic 
solids, whose properties are functions of the local state of deformation or stress. For example, the viscosity 
of a fluid may be a function of the local strain rates (velocity gradients), or the elastic modulus of a solid 
may be a function of the local state of stress or strain. The form of dependence in each case would be 
determined from experimental data. It is also possible for material properties, such as viscosity, thermal 
conductivity, or elastic modulus to be significantly dependent on temperature, and therefore not known 
in advance. Such problems cannot be solved by boundary element method which discretise only the 
boundary: some internal discretisation and analysis is required adequately to treat the non-linearity. At 


least some of the benefits of the boundary element approach may thereby be lost. 


7.5 Comparison with Other Methods 


Other methods for solving engineering continuum mechanics problems include finite element, finite 
difference and finite volume methods. They are what are often referred to as domain methods: the entire 
two- or three-dimensional domain of the problem has to be discretised, rather than just the boundary. 
The principal advantage of boundary element methods is the effective reduction in dimensionality of the 
numerical problem. At least for linear situations, a two-dimensional physical problem requires only a 
one-dimensional mesh of boundary line elements, while a three-dimensional physical problem requires a 
two-dimensional mesh of typically triangular or quadrilateral elements. Benefits of this simplification can 
be seen in the problems addressed in this book: meshes of line elements for two-dimensional problems 
are easily generated using only a small number of geometric parameters, leading to very concise data 


input files. Analyses of even quite complicated problems can be set up quickly and easily. 
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Although it has not been demonstrated in this book, boundary element meshes, particularly of quadratic 
elements, are not merely domain method meshes with the internal points removed. For a given level of 
boundary discretisation, boundary element methods normally give substantially greater accuracy. But, 
although far fewer algebraic equations have to be solved, they have to be treated as being fully populated 


with non-zero coefficients. The computational benefits, although still present, are less clearcut. 


Features of boundary element methods which are perhaps deterrents to their wider use by engineers 
are the mathematics involved. Not only are they arguably more complex than for domain methods, but 
also tend to be of unfamiliar types: for example, integral equations, singular fundamental solutions, and 
analyses often expressed in tensor notation. An aim of this book has been to try to simplify and clarify the 
mathematics involved, even when this is sometimes at the expense of more bulky and less mathematically 
elegant equations. A further aim has been to demonstrate that boundary element methods are worthy 


of serious consideration for engineering analyses, and are particularly straightforward to use. 
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Appendix A: Gaussian Quadrature 


The aim of this appendix is to explain the Gaussian quadrature method of numerical integration, 


particularly as applied to boundary elements. 


When a function cannot be integrated analytically, some form of numerical integration or quadrature 
must be used. Two of the simplest methods are the familiar trapezium and Simpson's rules. Given a 
function f(x), values of the function are obtained for a number of values of the independent variable x 
over the required range of integration, and an approximate value of the integral obtained from them. An 
important feature (and limitation) of these simple methods is that the values of the independent variable 
at which the function is sampled are often taken at constant intervals, and certainly little attention is 


paid to the optimum choice of sampling points. 


If there is no restriction on the positions at which the function can be evaluated, then there is a class 
of numerical integration procedures referred to as Gaussian quadrature, which is to be preferred. The 
advantage over the trapezium and Simpsons rules is that greater accuracy for a given number of function 


evaluations, and therefore of overall computational effort, can be obtained. 


{®@ 


=L é. 0 +1 


Figure A.1 Typical function to be integrated numerically 


Figure A.1 shows a typical function to be integrated: the area between the curve f(é) and the € axis 
between the limits € = —1 and € = +1 is required. It is convenient to transform the actual independent 
variable, such as x, into a dimensionless variable with these limits. This is the reason for introducing the 
local intrinsic co-ordinate § for quadratic boundary elements. Assuming that the function is evaluated 


at n points, known as the Gauss points, the required integral can be expressed approximately as 


[7 F@Oae ~ WE, Wy f Ey) (A.1) 
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where w, the are weighting factors. Simpson’ rule (defined in Equation 3.49) is actually a particular 
form of this result in which, for a single application of the rule 


§ =-1 §2 =0 3 =+1 (A.2) 


Given the freedom to choose the §,, however, the values of these and the corresponding weighting 
factors can be optimised for the integration of polynomial functions to give better accuracy. These are 
now well established and are available for a wide range of n, the number of Gauss points (for example, 
Stroud & Secrest 1966). 


For any numerical integration, choosing the number of Gauss points involves a compromise between 
accuracy and the amount of computation involved. In boundary element methods the number of points 
is normally at least n = 4. In this book a higher value of n = 8 is preferred, for which the values of € 


and W are 


€g = —&, = 0.9602898564 ww, = weg = 0.1012285362 
€7 = —& = 0.7966664774 W2 = W7 = 0.2223810344 (A.3) 
€ = —& = 0.5255324099 W3 = We = 0.3137066458 


& = —&, = 0.1834346424 wy, = ws = 0.3626837833 
to ten decimal places. 


Gaussian quadrature formulae have also been worked out for functions which involve particular forms of 
functions as multipliers, including cases of multipliers which are singular within the range of integration. 


One such case which is very relevant to boundary element methods is 


Jo In () fd ~ LP-1 wf (mp) (A.A) 
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the numerical values for n = 8 being 


ni = 0.013320244 
ni = 0.079750429 
n} = 0.197871029 


ni, = 0.354153994 


ni = 0.529458575 
ne = 0.701814530 
ni = 0.849379320 


ni = 0.953326450 


w; = 0.164416605 
W> = 0.237525610 
W3 = 0.226841984 


w4 = 0.175754079 


ws = 0.112924030 
we = 0.057872211 
w7 = 0.020979074 


Weg = 0.003686407 


(A.5) 


Reference 


Stroud, A.H. & Secrest, D. 1966, Gaussian Quadrature Formulas, Prentice-Hall. 
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Appendix B: Gaussian Elimination 


The aim of this appendix is to explain Gaussian elimination, which is a direct method for solving sets 


of simultaneous linear algebraic equations of the general form 


QA44X4 + a42X2 + 


a2{X41 + a22X2 + 


+ AtnXn = by 


+ aA2nXyn = b2 


(B.1) 
An1X1 + An2X2 + + AnnXn = Dy 
where X1,X2,"***,X, are the unknowns, and the coefficients a;; and b; are all known constants. This 
set of equations can be expressed in matrix form as follows 

ay. U2 Qin f~1 by 
Q21 22 A2n || *2 by 

ie ae = (B.2) 
Ant Qn2 Ann Xn by 

[A][x] = [5] (B.3) 


The unknowns are successively eliminated by algebraic manipulation. The first equation can be used 


to eliminate x, from the remaining n —1 equations. The modified second equation is then used to 


eliminate x2 from the remaining n — 2 equations, and so on until the last equation contains only x,. 
Thus may be found, followed by all the other unknowns, by back substitution. Let the coefficients of [A] 


and [b] shown in Equations B.1 and B.2 be given the notation ae and ae After the nth elimination, 


the modified coefficients are 


(k+1) _ 


z Uk), 
ij 


a ae — Oa,; 


(K+) _ p00 _ gp 
bh) — p® _ gp! 
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(B.4) 


where @ = an a, i=k+1,k+2,..,n andj =k,k + 1,...,n. Note that the column vector [b] 


is treated just like a column of [A], and advantage can be taken of this fact to simplify the computer 


programming of the process. The final set of equations is 


aX x + a x2 afr ees gedgeee te an x, = be 
2 2 2 
a x, Pianta ax, = p§ ) 


(B.5) 


a”) Xn, = pe 


Expressed in matrix terminology, the elimination process triangularises [A]. The unknowns are obtained 


in reverse order 
Xn = DM” sa? (B.6) 


n 


= (1-1 as} so? 
j=it1 


where! =n—1,n—-2,....,1. 


A great many arithmetic operations are involved in solving large sets of equations by elimination. Any 
errors introduced, such as roundoff errors caused by representing numbers with only a finite number 
of significant figures, tend to be magnified and may become unacceptably large. Equations B.4 show 
that the elimination process involves many multiplications by the factors @. In order to minimize the 
effects of roundoff errors, these factors should be made as small as possible, and certainly less than one. 
Therefore, the pivotal coefficient a“? should be the largest coefficient in the leading column of the 


kk 
remaining submatrix 


lat? | > lal?) fori=k+1,k+2,..,0 (B.7) 


This condition also helps to avoid division by zero in Equations B.6, and can be achieved by a process 
known as partial pivoting. Immediately before each elimination, the leading column is searched for the 
largest coefficient. By interchanging equations, this can be made the pivotal coefficient to satisfy Equation 
B.7. The idea of partial pivoting can be extended to searching the whole of the remaining submatrix 
for the largest coefficient. Such complete pivoting involves interchanging both rows and columns, and 
is more difficult to program. Since it offers only modest advantages in terms of accuracy over partial 
pivoting, it is rarely used. 
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Another refinement which helps to minimize the effects of roundoff errors is to scale the equations to 
make their coefficients similar in magnitude. One way to do this is to normalize each equation so that 
the largest coefficient in each row of [A] is of magnitude one. Scaling can be particularly important 


when corresponding coefficients in different equations differ by several orders of magnitude. 


Any given set of equations may provide either duplicate (for example, two equations identical) or 
inconsistent information, and cannot be solved. In either case, the matrix of coefficients [A] is said to 
be singular, and its determinant is zero. Such a condition can be detected during partial pivoting if it 
is impossible to find a non-zero pivotal value at some stage of the elimination process or at the start of 
the back substitution process. What is more difficult to detect, however, is when a set of equations is 
nearly singular or ill-conditioned. This arises when, for example, two equations provide very nearly the 
same information about the unknowns. Alternatively, the effect of roundoff errors may be to make what 
is a singular set of equations apparently ill-conditioned, in that rather than a zero pivotal coefficient, a 
very small value is obtained. Indeed, the detection of a very small pivotal coefficient can be used as a 
convenient test for either a singular or very ill-conditioned set of equations, but without being able to 
distinguish between them. By very small pivotal coefficient is meant a value very small in relation to the 
magnitudes of the coefficients of the matrix at the start of the elimination process, which are of order 


unity if scaling has been applied. 


The following subprogram ELIMIN provides a Fortran implementation of Gaussian elimination. 


SUBROUTIN 


r. 


ELIMIN (A, X,NEQN, NROW, NCOL, IFLAG) 


r. 


! SUBROUTINE FOR SOLVING SIMULTANEOUS LINEAR EQUATIONS BY GAUSSIAN 


! ELIMINATION WITH PARTIAL PIVOTING. 


REAL :: A(NROW,NCOL) , X (NROW) 


! INITIALIZE ILL-CONDITIONING FLAG. 


IFLAG=0 


! SCALE EACH EQUATION TO HAVE A MAXIMUM COEFFICIENT MAGNITUDE OF UNITY. 


JMAX=NEQN+1 


ach equation in turn: DO I=1,NEQN 


El 


AMAX=0. 


Search for maximum: DO J=1,NEQN 
ABSA=ABS (A(I,J) ) 
IF (ABSA > AMAX) AMAX=ABSA 


eal 


ND DO Search for maximum 
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Scale coefficients: DO J=1,JMAX 
A(1I,J)=A(1I,J) /AMAX 


END DO Scale coefficients 


END DO Each equation in turn 


a 


!  COMMENCE ELIMINATION PROCESS. 


Eliminate each variable in turn: DO K=1,NEQN-1 


! SEARCH LEADING COLUMN OF THE COEFFICIENT MATRIX FROM THE DIAGONAL 


! DOWNWARDS FOR THE LARGEST VALUE. 


r 


IMAX=K 


QN 
IF (ABS (A(I,K)) > ABS (A(IMAX,K))) IMAX=I 


1B 


Search for largest value: DO I=K+1,N 


eal 


ND DO Search for largest value 


! IF NECESSARY, INTERCHANGE EQUATIONS TO MAKE THE LARGEST COEFFICIENT 


! BECOME THE PIVOTAL COEFFICIENT. 


IF(IMAX /= K) THEN 
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Interchange coefficients: DO J=K, JMAX 
ATEMP=A (K, J) 

(K, J) =A (IMAX, J) 

(IMAX, J) =ATEMP 


END DO Interchange coefficients 


! ELIMINATE X(K) FROM EQUATIONS (K+1) TO NEON, 


Appendix B: Gaussian Elimination 


FIRST TESTING FOR 


| EXCESSIVELY SMALL PIVOTAL COEFFICIENT (ASSOCIATED WITH A SINGULAR 


! OR VERY ILL-CONDITION 


Gl 


D MATRIX). 


IF (ABS (A(K,K)) < 1.E-5) THEN 


IFLAG=1 
RETURN 


END IF 


Bach of remaining equations: DO I=K+1,N 
FACT=A (I,K) /A(K,K) 

Modify coefficients: DO J=K, JMAX 
A(I,J)=A(1I,J) -FACT*A(K,J) 

END DO Modify coefficients 


END DO Bach of remaining equations 


eal 


ND DO Eliminate each variable in turn 


eal 
10 
Zz 


! SOLVE THE EQUATIONS BY BACK SUBSTITUTION, FIRST TESTING 


! FOR AN EXCESSIVELY SMALL LAST DIAGONAL COEFFICIENT. 


IF (ABS (A(NEQN, NEQN)) < 1.E-5) THEN 


IFLAG=1 
RETURN 


a 
z 


D IF 


X (NEQN) =A (NEON, JMAX) /A (NEON, NEQN) 


Then each unknown in turn backwards: DO I=NEQN-1,1,-1 


UM=A (I, JMAX) 


S 
Sum products: DO J=I+1,NEQN 
SUM=SUM-A (I,J) *X (J) 

END DO Sum products 

X (I) =SUM/A(I,T) 


END DO Then each unknown in turn backwards 


7) 
irl 


ETURN 


END SUBROUTINE ELIMIN 
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The arguments of ELIMIN include the array of equation coefficients, A, which incorporates the vector 
of constants [b] as its last column, and the solution vector, X. The argument NEQN enters the number 
of equations to be solved, while IFLAG returns an integer flag value to the calling program to warn of 


a singular or very ill-conditioned coefficient matrix. 


The first action of the program is to scale the coefficients of the equations to give a maximum coefficient 
magnitude of unity in each row of [A]. The elimination process is then started: each equation, with 
the exception of the last, is used in turn to eliminate the corresponding unknown from the subsequent 
equations. The current elimination number is given by K, and is equivalent to k in Equations B.4. Before 
performing the necessary eliminations with a particular equation, however, a search is made down the 
leading column of the remaining submatrix to find the coefficient of greatest magnitude, as defined 
by Equation B.7. The search technique locates the row number of the largest coefficient, IMAX, by 
first assuming that it corresponds to the coefficient on the diagonal, and only changing this if a larger 
coefficient is found. If the pivotal coefficient is not the largest, the relevant equations are interchanged 
by interchanging all their coefficients. In computing terms, this movement of data between storage 
locations is inefficient. While it could be avoided by keeping a record of the revised order in which the 
equations are to be considered, this makes the program significantly more difficult to understand. For 


present purposes, some computational efficiency is sacrificed in the interests of clarity. 


Despite the search for the largest coefficient, it is still possible for the pivotal coefficient to be extremely 
small or zero if the coefficient matrix is singular or very ill-conditioned. Bearing in mind that the equations 
were scaled initially, an appropriate test of magnitude is |a;,;,| < 10~>. If this condition is satisfied at any 
stage of the elimination process, the problem is rejected. Rejection is indicated by setting the value of 
IFLAG to one, which may be detected by the calling program. Ifthe partial pivoting is successful, however, 
the eliminations defined by Equations B.4 are performed, with the variable FACT being used to store 
the values of the factor @. After testing the magnitude of the last diagonal coefficient (a in Equations 


B.5), the back substitutions defined in Equations B.6 are performed to find the required solutions. 
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Appendix E: Matlab Version of 
Quadratic Boundary Element 
Program for Plane Elastic Problems 


The aim of this appendix is to present a Matlab version of the program BEM2EQ which is described in 


detail in its Fortran form in Chapter 6. Matlab is now very widely used by engineers for many purposes. 


The Matlab version is as far as possible a literal translation of the Fortran. This means that the detailed 
explanations provided in Chapter 6 are also applicable to the Matlab script. On the other hand, no 
attempt is made to take advantage of special built-in features of Matlab such as simplified handling of 


matrices and vectors, and the solution of sets of linear equations. 


Fortran is a compiler which converts the source code into machine code before computation starts. In 
contrast, Matlab is an interpreter which goes through the code line-by-line, translating and executing 
each line as it goes. Strictly speaking therefore, Matlab works with scripts rather than programs, although 
scripts are often loosely referred to as programs. A consequence of this difference between a compiler and 
an interpreter is that for large computationally intensive problems Fortran generally executes noticeably 


more quickly than Matlab. 
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The same function and variable names are used in the Fortran and Matlab programs: upper case names 
in Fortran become lower case in Matlab. The definitions of the names given in Some Program Variable 
Names at the beginnings of both Parts of the book are equally valid for both. DO loops in Fortran become 
“for” loops in Matlab, and functions are called in somewhat different ways. Also, STOP in Fortran is 
replaced by an error message to the computer screen in Matlab. The main difference, however, is to be 
seen in the details of the handling of input and output, although the same approach of inputting data 
from one pre-defined file and outputting to other files is adopted. In particular, the required input DATA 


files used in Fortran are unchanged, and the output files are effectively identical. 


The Matlab version of BEM2EQ, stored as the file BEM2EQ.m, is as follows. 


function BEM2EQ (varargin) 


Q 
© 


ae 


PROGRAM FOR SOLVING TWO DIMENSIONAL ELASTIC STRESS ANALYSIS PROBLEMS 


oe 


BY THE BOUNDARY ELEMENT METHOD USING QUADRATIC ELEMENTS. 


Q 
© 


clear global; clear functions; 


global fid5 fide fid7 


shareddata2eq 
fid5=fopen('DATA', 'r'); 
fid6é=fopen('RESULTS', 'wt'); 


fid7=fopen ('MESHRES', 'wt'); 


3 
io) 


% DEFINE THE MAXIMUM PROBLEM SIZE 


U 


ERMITTED BY THE ARRAY DIMENSIONS. 


maxnel=250; 
maxnnp=500; 
maxnb=10; 
maxnpc=20; 


Q 
© 


r. 


% INPUT THE PROBLEM TITLE AND TYPE, ALSO MATERIAL PROPERTIES. 


intitle; 


Q 
© 


% INPUT AND GENERATE THE MESH DATA. 


meshq; 


eo 
io) 


% OUTPUT THE 


< 


ESH DATA. 


mshout; 


3 
io) 
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% EVALUATE AND STORE VALUES OF THE SHAPE FUNCTIONS AND THEIR DERIVATIV. 


Gl 
wn 


% AT THE GAUSS POINTS AND NODES. 


zr 


shape; 


% 


% INPUT, PROCESS AND OUTPUT THE BOUNDARY CONDITIONS. 


bes; 


% 


[I 


% FORM THE COEFFICIENT MATRIX AND APPLY TH 


BOUNDARY CONDITIONS. 


frmtrx; 


% 


% SOLVE THE LINEAR EQUATIONS. 
neqn=2*nnp; 
[uv, iflag]=elimin(a,negqn) ; 
if (iflag == 1) 
fprintf (fid6, ['\n', 


"MATRIX ILL-CONDITIONING DETECTED IN EQUATION SOLVER','\n']); 


error ('MATRIX ILL-CONDITIONING DETECTED IN EQUATION SOLVER") ; 


end; 


% 


G 


% OUTPUT THE NODAL DISPLAC 


es 


(MENTS, ELEMENT STRESSES AND FORCES ON THE 
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% BOUNDARY SEGMENTS. 


output; 


% 


end Sprogram BEM2PQ 


function intitle (varargin) 


foe) 


oP 


SUBPROGRAM TO INPUT PROBLEM TITLE AND WHETHER PLANE 


= 


% STRESS OR PLANE STRAIN. ALSO THE MATERIAL PROPERTIES. 


global nu e fid5 fid6é 


% 


% INPUT THE PROBLEM TITLE. 


title=fgetl (fid5) ; 


r 


fprintf (fid6é, ['QUADRATIC BOUNDARY EMENT SOLUTION FOR', 


' TWO DIMENSIONAL ELASTIC PROBLEM','\n','\n','%Ss','\n'], title); 


% INPUT THE PROBLEM CASE TYPE: 


% "STRESS' DEFINES PLANE STRESS 
% 'STRAIN' DEFINES PLANE STRAIN 


fo) 


THE DEFAULT IS PLANE STRESS. 


pcase=fgetl (fid5) ; 


fprintf (fid6é, ['\n', 'UNDER PLANE ','%Ss',' CONDITIONS','\n'],pcase) ; 


% 


& INPUT AND OUTPUT YOUNG'S MODULUS AND POISSON'S RATIO. 
line=fgetl (fid5) ; 
vec=str2num(line); 


e=vec(1); nu=vec(2); 


fprintf (fid6é, ['\n', 'YOUNGS MODULUS = ','%12.4e', 
’ POISSONS RATIO = ','%5.3f£','\n'],e,nu); 


% 


% MODIFY PROPERTIES IF CASE IS ONE OF PLANE STRAIN. 


strain='STRAIN'; 


if (strcmpi(pcase,strain) == 1) 
e=e/(1.-nu%*2); 
nu=nu/ (1.-nu) ; 
end; 
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return; 


end function intitle 


function meshq(varargin) 


% 


% SUBPROGRAM TO READ IN AND GENERATE THE GEOMETRIC DATA FOR A MESH OF 


% QUADRATIC ELEMENTS. 


% 


global node maxl ynode xnode nnp nel neend m3 ml angstore yeend xeend 


isegelem mfirst isegend ysend xsend mlast maxnel nelb nsegb 


nsegtot nbound maxnb fid5 fid6 


% 


% INPUT THE NUMBER OF SEPARATE BOUNDARIES. 


line=fgetl (fid5) ; 
nbound=str2num (line) ; 


% 


$ TEST THE NUMBER OF BOUNDARIES. 
if (nbound < 1 || nbound > maxnb) 


fprintf (fidé, ['\n', 'NBOUND =','34i', 


'  OUTSID 


I 


E 


PERMITTED RANGE 1 TO', '%4i','\n'],nbound,maxnb) ; 
error('nbound error'); 


end; 


% 


% FOR EACH BOUNDARY IN TURN INPUT THE NUMBER OF SEGMENTS. 


nel=0; 
ieend=0; 
nsegtot=0; 
mmaxb=0; 
for ibound=1:nbound; % each boundary in turn 
nelb (ibound) =0; 
mminb=mmaxb+1; 
line=fgetl (fid5) ; 
nsegb (ibound) =str2num (line) ; 


nsegtot=fix (nsegtottnsegb (ibound) ); 


foe) 


6 TEST THE NUMBER OF SEGMENTS. 


if(nsegtot < 1 || nsegtot > maxnel) 


fprintf (fid6é, ['\n', 'NS 


Download free eBGUKS at bookboow dam ED RANGE 1 TO','%S6i','\n'],nsegtot,maxnel) ; 


Gl 


GTOT =', 'S6i', 
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error('nsegtot error'); 


end; 


% INPUT THE CARTESIAN GLOBAL COORDINATES OF THE END POINTS OF THI 


Gl 


% SEGMENTS. TAKE THE END POINTS CONSECUTIVELY, K 


zr 


‘EPING THE DOMAIN 


6 TO THE LEFT OF THE DIRECTION OF NUMBERING. 


% Each line of input data is assumed to contain an arbitrary number 
% of pairs of x,y coordinates. 
isend=0; 
while isend < nsegb(ibound) ; 
line=fgetl (fid5) ; 
vec=str2num(line) ; 
veclength=length (vec) ; 
np=veclength/2; 
for ip=l:np; 
isend=isendt+1; 


xsend(isend) =vec(2*ip-1); 


ysend(isend) =vec(2*ip); 


end; 
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oP 


DEFINE THE FIRST END POINT ON TH 


r 


CURRENT BOUNDARY. 


jeend=ieend+1; 


xeend(ieend) =xsend (1); 
yeend (ieend) =ysend (1); 


% FOR EACH OF THE SEGMENTS (BETWEEN ENDS 1 AND 2, 2 AND 3, ETC.) 


% INPUT TH 


ira 


RADIUS OF CURVATURE (+VE FOR CONVEX WITH CENTRE OF 


% CURVATUR 


Gl 


INSIDE DOMAIN, -VE FOR CONCAVE), THE NUMBE 


Pe) 
oO 
Hy 


% ELEMENTS IN THE SEGMENT, AND THE LENGTH RATIO BETWEEN SUCCESSIVE 


% ELEMENTS IN THE DIRECTION OF NUMBERING. 


isegmax=nsegtot; 


isegmin=isegmax-nsegb (ibound) +1; 


for iseg=isegmin:isegmax; % each segment in turn 
line=fgetl (fid5) ; 
vec=str2num(line) ; 


rseg=vec(1); nelseg=vec(2); ratseg=vec (3); 


% 


% FIND AND TEST THE NUMBER OF ELEMENTS SO FAR. 


nel=fix (nelt+tnelseg) ; 


nelb (ibound) =fix (nelb (ibound) +nelseg) ; 


if(nmel < 1 || nel > maxnel) 


fprintf (fid6, ['\n', 'NEL =','%6i',' OUTSIDE PERMITTED! 


' RANGE 1 TO','%6i','\n'],nel,maxnel); 
error('nel error'); 


end; 


oe 


foe) 


FIRST AND LAST ELEMENTS ON CURRENT SEGMENT. 
mlast (iseg) =fix (nel) ; 


mfirst (iseg) =fix (nel-nelsegt1); 


mmaxb=mmaxb+nelseg; 


oe 


oe 


COORDINATES OF THE FIRST END POINT OF THE SEGMENT. 


isend=iseg-isegmint1; 


xfirst=xsend(isend) ; 


yfirst=ysend (isend) ; 


foe) 


ol? 


COORDINATES OF THE LAST END POINT OF THE 


= 
n 


EGMENT . 


isend=isend+1; 
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if(iseg == isegmax) 
isend=1; 

end; 

xlast=xsend(isend) ; 


ylast=ysend (isend) ; 


oe 
@ 
Ss 
fe 
= 
s 


EMENT DATA FOR A STRAIGHT SEGMENT. 


fo) 


% DEFINE THE ELEMENT END POINT COORDINATES ON THE SEGMENT. 


for m=l:nelseg; % each element in turn 


ieend=ieend+1; 


isegend (ieend) =fix (iseg) ; 
if(ratseg == 1.) 
xeend (ieend) =xfirst+ (xlast-xfirst) * (m) / (nelseg) ; 
yeend (ieend) =yfirst+ (ylast-yfirst) * (m) / (nelseg) ; 
end; 
if(ratseg ~= 1.) 
xeend (ieend) =xfirst+ (xlast-xfirst) * (1.-ratseg*m) / 
(1.-ratseg*nelseg) ; 
yeend (ieend) =yfirst+ (ylast-yfirst) * (1.-ratseg*m) / 
(1l.-ratseg*nelseg) ; 
end; 


fo) 


end; % each element in turn 


6 DEFINE THE ELEMENT NODES AND COORDINAT 


GI 
nN 


ijeend=ieend-nelseg; 


for ielseg=l:nelseg; % each element in turn 
ieend=ieend+1; 
m=mfirst (iseg)+tielseg-1; 
i11=2*m-1; 
i2=il1tl; 
13=114+2; 


if(iseg == isegmax && ielseg == nelseg) 


i3=node(mminb,1); 
end; 
node (m, 1) =fix (il); 


node (m, 2) =fix (12); 


node (m, 3) =fix (13); 


isegelem (m) =fix (iseg) ; 
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if(iseg == isegmin && ielseg == 1) 
xnode (il) =xeend(ieend-1); 
ynode (i1) =yeend (ieend-1) ; 


end; 


xnode (i3)=xeend(ieend) ; 


ynode (13) =yeend (ieend) ; 
xnode (12) =0.5* (xnode (il) +xnode(i3)); 


ynode (12) =0.5* (ynode (il) +ynode(i3)); 


foe) 


STORE THE NUMBERS OF THE ADJACENT ELEMENTS. 


oe 


ml (m) =fix (m-1) ; 
m3 (m) =fix (m+1) ; 
end; % each element in turn 


end; 


% GENERATE ELEMENT DATA FOR A SEGMENT IN THE FORM OF A CIRCULAR ARC. 


if(rseg ~= 0.) 


fo) 


LOCATE THE CENTRE OF THE ARC. 


ol? 


xmid=(xfirsttxlast) /2.; 


> Apply now 


REDEFINE YOUR FUTURE 
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PROGRAM 2015 
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ymid=(yfirstt+tylast)/2.; 


alseg=sqrt ((xlast-xfirst) *2+(ylast-yfirst) *2); 


alperp2=rseg*2-(alseg/2.)%*2; 


if (abs (alperp2) 


alperp2=0.; 


end; 


< 1.0e-6*rseg%2) 


if(alperp2 < -1.0e-6*rseg%2) 


fprintf (fid6, ['\n', 'DATA ERROR FOR SEGMENT NUMBER','36i', 


"\n', 'NOT POSSIBLE TO CREATE A CIRCULAR ARC','\n'],iseg); 


error('circula 


end; 


alperp=sqrt (alpe 


uvflx=(xlast-xfirs 


uvfly=(ylast-yfirs 


fact=1.; 


if(rseg < 0.) 


r arc error'); 


rp2); 
t) /alseg; 


t) /alseg; 


ENDED THERE BY THE SEGMENT. 


fact=-1.; 
end; 
xcent=xmid-alperp*uvfly* fact; 
ycent=ymidtalperp*uvilx* fact; 
% FIND THE ANGLE SUBT 
if(alperp ~= 0.) 


angseg=2.*atan(alseg*0.5/alperp) ; 


end; 


if(alperp == 


angseg=pi; 


end; 


fo) 


oe 
iw) 


EFINE 


THE 


G 


EL 


EM 


ENT 


angfir=atan2 (yfirs 


t-ycent, xfirst-xcent) ; 


angstore (ieend) =0.; 


for m=l1:nelseg; 


jeend=ieend+1; 


isegend(i 


nd) 


% each element in turn 


if(ratseg == 


en 


=fix (iseg) ; 


:) 


ang=angseg* (m) / (nelseg) ; 


d; 


if(ratseg ~= 1.) 


END POINT COORDINATES ON THE 


EGMENT. 


ang=angseg* (1.-ratseg*m) /(1.-ratseg*nelseg) ; 
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end; 
if(rseg < 0.) 
ang=-ang; 


end; 


xeend (ieend) =xcenttabs (rseg) *cos (angfir+ang) ; 


yeend (ieend) =ycentt+abs (rseg) *sin(angfir+ang) ; 
angstore (ieend) =ang; 


fo) 


end; % each element in turn 


%& DEFINE THE ELEMENT NODES AND COORDINAT 


Gl 
Nn 


ieend=ieend-nelseg; 


for ielseg=l:nelseg; % each element in turn 
ieend=ieend+1; 
m=mfirst (iseg)+tielseg-1; 
i11=2*m-1; 
i2=il1tl; 


13=1i1+2; 


if(iseg == isegmax && ielseg == nelseg) 
i3=node(mminb,1); 

end; 

node (m, 1) =fix (i1); 

node (m, 2) =fix (12) ; 


node (m, 3) =fix (13); 


isegelem(m) =fix (iseg) ; 

if(iseg == isegmin && ielseg == 1) 
xnode (il) =xeend(ieend-1) ; 
ynode (i1)=yeend (ieend-1) ; 


end; 


xnode (i3)=xeend(ieend) ; 


ynode (13) =yeend(ieend) ; 


ang=0.5* (angstore (ieend-1) +angstore (ieend) ); 
xnode (i2) =xcent+abs (rseg) *cos (angfirtang) ; 


ynode (i2) =ycent+abs (rseg) *sin(angfir+ang) ; 


% STORE THE NUMBERS OF THE ADJACENT ELEMENTS. 


a9 


ml (m) =fix (m-1) ; 
m3 (m) =fix (m+1) ; 


end; % each element in turn 


end; 
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end; % each segment in turn 


%& ADJACENT ELEMENTS FOR END ELEMENTS OF THE BOUNDARY. 


ml (mminb) =fix (mmaxb) ; 


ix (mminb) ; 


m3 (mmaxb) = 
end; % each boundary in turn 
neend=fix (ieend) ; 


nnp=fix (nel*2) ; 


% DETERMINE THE MAXIMUM DIMENSION OF THE SOLUTION DOMAIN. 


maxl=0.; 


for i=l:nnp; % each node in turn 
for j=l:nnp; % each other node in turn 
dist=sqrt ( (xnode (i) -xnode (j) ) *2+ (ynode (i) ~ynode(j))*2); 
if (dist > maxl) 
maxl=dist; 
end; 
end; % each other node in turn 


end; % each node in turn 
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return; 


end Sfunction meshq 


function mshout (varargin) 


% 


\ 


% SUBPROGRAM TO WRITE OUT THE MESH DATA. 


= 


% 


global node maxl ynode xnode nnp nel fid7 


% 


% OUTPUT TH 


r 


NUMBERS OF ELEMENTS AND NODES, ALSO THE NODAL 


% COORDINATES. 


fprintf (fid7, ['GEOMETRIC DATA FOR THE MESH','\n','\n', 


' NUMBER OF ELEMENTS =','36i','\n','\n', 


: NUMBER OF NODAL POINTS =','%S6i','\n','\n', 


"COORDINATES OF THE NODAL POINTS','\n','\n', 


; I X Y , 
' I xX Y '],nel,nnp); 
for i=l:nnp; 
if (mod(i,2) == 1); 
fprintf (fid7, ['\n','S6i','S12.4e','S12.4e'], 
i,xnode(i),ynode(i)); 
end; 
if (mod(i,2) == 0); 
fprintf (fid7, ['S61i','%S12.4e','%12.4e'], 


i,xnode(i),ynode(i)); 


% OUTPUT TH 


29 


ELEMENT NODE NUMBERS. 


fprintf (fid7, ['\n','\n', 'ELEMENT NODE NUMBERS','\n','\n', 
' M ND1 ND2 ND3', 


: M ND1 ND2 ND3']); 


for m=l1:nel; 


if(mod(m,2) == 1) 
fprintf (fid7,['\n',' ","S5i'],m); 
end; 
if(mod(m,2) == 0) 
fprintf (fid7, [' ","S5i'],m); 
end; 
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for in=1:3; 
forintf (fid7,'S5i',node(m,in)); 


end; 


% SCALE THE NODAL POINT COORDINATI 


ea 
n 


for i=l:nnp; % each node in turn 
xnode (i) =xnode (i) /maxl; 
ynode (i) =ynode (i) /maxl1; 

end; % each node in turn 

return; 


end Sfunction mshout 


function bcs (varargin) 


% 


% SUBPROGRAM TO INPUT, PROCESS AND OUTPUT TH 


Gl 


BOUNDARY CONDITIONS. 


% 


global nodepe idirpe nnp nbcpe nbcu store node tmy ty tmx tx 


ibcn ibce maxl velem v uelem u m3 ml nel mlast mfirst sigsnseg 
signnseg sigsn signn unmx unmy nsegtot isegbc nbcs vseg useg 


maxnpc maxnel nbct ungy ungx fid5S fid6 


% 


ae 


zi 
t= 
4 


% FIRST FIND THE UNIT NORMAL COMPON 


Gl 


NTS AT THE 


EMENT NOD! 


Gl 
wn 


for m=l:nel; % each element in turn 


for in=1:3; % each element node in turn 
igauss=in+10; 
it=1; 
ic=1; 
jacobi(m,igauss,it,ic); 
unmx (m, in) =ungx; 
unmy (m, in) =ungy; 
end; % each element node in turn 


(a) 


end; % each element in turn 


% 


fo) 


INPUT THE NUMBERS OF SEGMENTS SUBJECT TO EACH TYP 


Gl 


OF BOUNDARY 


oe 


CONDITION. ALSO THE NUMBER OF POINT CONSTRAINTS. 


oe 


NBCU - PRESCRIBED DISPLACEMENTS. 


fo) 


NBCS - NON-ZERO PRESCRIBED STRESSES. 
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% ANY SEGMENT NOT INCLUDED IS ASSUMED TO BE SUBJECT TO ZERO 

% STRESSES. 

% NBCPC - NODAL POINT DISPLACEMENT CONSTRAINTS. 

line=fgetl (fid5) ; 

vec=str2num(line); 

nbcu=vec(1); nbcs=vec(2); nbcpc=vec (3); 

% TEST THESE BOUNDARY CONDITION NUMBERS. 

nbct=fix (nbcutnbcs) ; 

if(nbcu < 0 || nbcu > maxnel || nbcs < 0 || nbcs > maxnel 

|| nbct < 0 || nbct > maxnel) 
fprintf (fid6é, ['\n', 'NBCU =','S6i',' NBCS =','%S6i','\n', 
: NBCT =!,'%S6i',' OUTSIDE PERMITTED RANGE 0 TO', 


'S6i','\n'],nbcu,nbcs,nbct,maxnel); 
error('numbers of boundary conditions error'); 
end; 
if(nbcpe < 0 || nbcpce > maxnpc) 


fprintf (fid6é, ['\n','NBCPC =','S4i', 


' OUTSIDE PERMITTED RANGE 0 TO','%S3i','\n'],nbcpc,maxnpc) ; 
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error('number of point constraints error'); 


end; 


% 


\ 
= 


6  INITIALIS 


r 


THE BOUNDARY CONDITION STORAGE ARRAYS. 


for m=l:nel; % each element in turn 


ibce (m)=2; 


for in=1:3; % each element node in turn 
signn(m,in)=0.; 
sigsn(m,in)=0.; 
tmx (m,in)=0.; 
tmy (m,in)=0.; 
end; % each element node in turn 


fe) 


end; % each element in turn 


% 


% INPUT, STORE AND OUTPUT THE PRESCRIBED DISPLACEMENT CONDITIONS. 


if (nbcu > 0) 
for ibcu=1:nbcu; 
line=fgetl (fid5) ; 
vec=str2num(line) ; 
isegbc (ibcu)=vec(1); useg(ibcu) =vec(2); vseg(ibcu)=vec (3); 


end; 


fprintf (fid6é, ['\n', 'PRESCRIBED DISPLACEMENT BOUNDARY CONDITIONS','\n']); 
for ibcu=l:nbcu; % each segment with prescribed displacements 
iseg=isegbc (ibcu) ; 


if(iseg < 1 || iseg > nsegtot) 


fprintf (fid6é, ['\n', 'ISEG = ','%6i',' OUTSIDE PERMITTED RANGE' 
' 1 T0','S6i','\n'],iseg,nsegtot) ; 

error('iseg error'); 

end; 

for m=mfirst (iseg) :mlast(iseg); % each element on current segment 
ibce (m)=1; 
uelem(m)=useg (ibcu) ; 
velem(m)=vseg (ibcu) ; 


end; % each element on current segment 


oe 


fprintf (fid6é, ['\n','U =', '$12.4e',' V =','%12.4e', 


' ON ELEMENTS','34i',' TO’, 34a", \n"' Jy 


useg (ibcu) , vseg(ibcu) ,mfirst (iseg),mlast (iseg) ); 
end; % each segment with prescribed displacements; 


end; 
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% INPUT, STORE AND OUTPUT THE PRESCRIBED STRESS CONDITIONS. 


if(nbcs > 0) 
for ibcs=l1:nbcs; 
line=fgetl (fid5) ; 
vec=str2num (line) ; 
isegbc (ibcs)=vec(1); signnseg(ibcs)=vec(2); sigsnseg(ibcs) =vec (3); 


end; 


fprintf (fide, ['\n','PRESCRIBED STRESS BOUNDARY CONDITIONS','\n']); 


for ibcs=l:nbcs; % each segment with prescribed stresses 
iseg=isegbc(ibcs); 


if(iseg < 1 || iseg > nsegtot) 


fprintf (fidé, ['\n','ISEG = ','%6i',' OUTSIDE 


U 


ERMITTED RANGE' 


1 TO', 'S6i','\n'],iseg,nsegtot) ; 
error('iseg error'); 
end; 
for m=mfirst (iseg) :mlast(iseg); %© each element on current segment 


ibce(m)=2; 


oe 


ol? 


FIND THE TRACTIONS AT THE ELEMENT NODES FROM THE PRESCRIBED STRESSES. 


foe) 


ALSO STOR 


r 


THE STRESSES AT THE ELEMENT NODES. 


for in=1:3; % each element node in turn 
tmx (m, in) =signnseg (ibcs) *unmx (m, in) -sigsnseg(ibcs) *unmy(m, in) ; 
tmy (m, in) =signnseg(ibcs) *unmy (m, in) +Sigsnseg(ibcs) *unmx(m, in) ; 
signn(m,in)=signnseg(ibcs) ; 
sigsn(m,in)=sigsnseg(ibcs) ; 


end; % each element node in turn 


end; % each element on current segment 


oe 


fprintf (fid6é, ['\n', 'SIGNN =','$12.4e',' SIGSN =','%12.4e', 
' ON ELEMENTS', '34i',' TO', '$4i','"\n'], 


signnseg(ibcs),sigsnseg(ibcs) ,mfirst (iseg) ,mlast (iseg)); 


end; % each segment with prescribed stresses 


end; 


eo 
) 


% ASSEMBLE THE VECTOR OF KNOWN VARIABLES, APPLYING SCALING. 


& FIRST INITIALISE THE NODAL DISPLACEMENTS AND TRACTIONS. 


for i=l:nnp; % each node in turn 


u(i)=0.; 
v(i)=0.; 
tx(i)=0.; 
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end; % each node in turn 


for m=l:nel; % each element in turn 


for in=1:3; % each element node in turn 
i=node(m,in); 


ibcn (i) =fix (ibce(m) ); 


% CHECK THE CONDITION APPLIED TO THE ADJACENT ELEMENT FOR AN ELEMENT 


% end NODE. IF THE CONDITION IS PRESCRIBED DISPLACEMENTS BUT THE 


% EXISTING CONDITION AT THE NODE IS NOT, THEN IMPOSE PRESCRIBED 


& DISPLACEMENTS. 

madj=m; 

if(in == 1) 
madj=m1 (m) ; 

end; 

if(in == 3) 
madj=m3 (m) ; 

end; 

if (ibce(madj) == 1 && ibcn(i) == 2) 
iben(i)=1; 


end; 
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% STORE KNOWN VARIABLES. 


if(v(i) == 0.) 


end; 
if(ibce(m) == 2 && ibcn(i) 
if(tx(i) == 0.) 
tx(i)=tmx(m,in)/e; 


end; 


end; 


end; 


end; % each element node in turn 


end; % each element in turn 


% 


% SCALE YOUNG'S MODULUS. 


estore=e; 
e=1.; 


% 
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a 


% INPUT, STORE AND OUTPUT NODAL POINT DISPLACEMENT CONSTRAINTS. 


\ 


6 IDIRPC STORES DIRECTIONS OF THE CONSTRAINTS - 1 FOR X, 2 FOR Y. 


% SUCH CONSTRAINTS SHOULD NOT REQUIRE POINT FORCES TO MAINTAIN THEM. 


if (nbcu == 0 && nbcpe < 3) 
forintf (fid6é, ['\n', 'INSUFFICI 


EA 


NT DISPLACEMENT BOUNDARY CONDITIONS', 


' OR POINT CONSTRAINTS','\n']); 


error('displacement constraint error'); 
end; 
if (nbcpc > 0) 
for ibcpc=1:nbcpc; 
line=fgetl (fid5) ; 
vec=str2num(line) ; 
nodepc (ibcpc) =vec(1); idirpe (ibcpc)=vec (2); 
end; 
for ibcpc=l:nbcpc; % each point constraint in turn 


if (nodepc(ibcpc) < 0 || nodepc(ibcpc) > nnp) 


fprintf (fid6, ['\n', 'POINT CONSTRAINT NODE','34i', 


' OUTSIDE RANGE 0 TO ','%4i','\n'],nodepc(ibcpc),nnp) ; 


error('point constraint node error'); 
end; 


if(idirpce(ibcpc) ~= 1 && idirpc(ibcpc) ~= 2) 


fprintf (fid6, ['\n', 'POINT CONSTRAINT DIRECTION', '%3i', 


' FOR NODE','34i',' OUTSIDE RANGE 1 TO 2','\n'], 


idirpce (ibcpc) ,nodepe (ibcpc) ); 
error('point constraint direction error'); 
end; 


if(idirpce(ibcpc) == 1) 


fprintf (fid6é, ['\n', 'NODE', 'S4i',' CONSTRAINED WITH U = 0O', 
"\n'],nodepc (ibcpc) ); 
end; 
if(idirpc(ibcpc) == 2) 
fprintf (fid6é, ['\n', 'NODE', 'S4i',' CONSTRAINED WITH V = 0', 


"\n'],nodepc(ibcpc) ); 
end; 
end; © each point constraint in turn 
end; 
return; 


end function bcs 
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function frmtrx(varargin) 


% 


oe 


SUBPROGRAM TO FORM THE COEFFICIENT MATRIX AND RIGHT HAND SIDE 


< 


ECTOR, 


oe 


MODIFIED TO SUIT THE BOUNDARY CONDITIONS. 


2 


global a idirpc nodepc nbcpc node v arowy u arowx ibcn estore tmy 


browy tmx browx ibce nel nnp 


% 


% DEFINE THE NUMBER OF COLUMNS IN TH 


zr 


EXTENDED COEFFICIENT MATRIX A. 


jJmax=2*nnptl1; 


% 


% FORM THE COEFFICIENT MATRIX A, AND THE RIGHT HAND SIDE 


< 


ECTOR B*T. 


for i=l:nnp; % take each node in turn as P 


% 


% INITIALISE THE TWO ROWS OF MATRIX A CORRESPONDING TO THE NODE. 


fo) 


for j=l:jmax; % each pair of coefficients in turn 


DO YOU WANT TO KNOW: 


What your staff really want? 


The top issues troubling them? 


How to make staff assessments 
work for you & them, painlessly? 


How to retain your 
top staff 


FIND OUT NOW FOR FREE 


Get your free trial 


Because happy staff get more done 


Download free eBooks at bookboon.com 


161 Click on the ad to read more 


Boundary Element Methods for Engineers: 
Part 2: Plane Elastic Problems 


fo) 
(o) 


EL 


\ 


G 


THE 


EM 


INITIALIS 


a 


% 


A AND B MATRICES. 


% each 


for m=l1:nel; 


for in=1:3; $% 
arowx (m 
arowx (m 
arowy (m 
rowy (m 
rowx (m 


rowx (m 


a 
b 
b 
b 


rowy (m 


browy (m 


% 


end; 


% each element in turn 


end; 


ENT CONTRIBUTIONS TO TH 


Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


ENT ROWS OF TH 


Gl 


CURR 


lement in turn 


each element node in turn 


each element node in turn 


LOOP TO INT 


SET UP THE EGRATE OV 


r 


ER EL 


EM 


KACH ENT IN TURN. 


% each element in t 


for m=l:nel; 


INT EL 


G 


EGRATE 


G 


THE 


ERN 


intker (i,m); 


% each element in turn 


end; 


oP 


foe) 


ice) 


EFFICI 


EVALUATE 


TH ENTS AT TH 


aiixx=0.; 
aiixy=0. 
aiiyx=0.; 


, 


aiiyy=0.; 


, 


ach 


% 


for m=l1:nel; 
for in=1:3; % 
aiixx=aiixx-arowx(m,2*in-1); 
aiixy=aiixy-arowx(m,2*in); 
ailiyx=aiiyx-arowy(m,2*in-1); 
aliyy=aiiyy-arowy(m,2*in); 

end; % 

each element in turn 


) 


a(2*i-1,2*i-1)=aiixx; 


% 


end; 


if (iben (i) 


a(2*i-1,2*1i)=aiixy; 
a(2*i,2*i-1)=aiiyx; 


a(2*i,2*i)=aiiyy; 
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EM 


CURR ENT. 


a 
Bi 


DIAGONAL OF MATRIX A. 
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each element node in turn 


162 


Boundary Element Methods for Engineers: Appendix E: Matlab Version of Quadratic Boundary 
Part 2: Plane Elastic Problems Element Program for Plane Elastic Problems 


% 


\ 


6 INITIALIS 


zr 


THE B*T VECTOR COEFFICIENTS. 


btx=0.; 
bty=0.; 
if(ibcn(i) == 1) 


btx=-aiixx*u(i)-aiixy*v(i); 


bty=-aiiyx*u(i)-aiiyy*v(i); 


end; 


foe) 


% APPLY THE BOUNDARY CONDITIONS TO THE CURRENT ROWS OF A AND B, BY 


ol? 


CONSIDERING EACH ELEMENT IN TURN. 


for m=l:nel; % each element in turn 


for in=1:3; % each element node in turn 


j=node (m,in); 


fo) 


IF DISPLACEMENTS ARE PRESCRIBED OVER THE ELEMENT, 


ol? 


INTERCHANGE THE A AND B COEFFICIENTS. 


if(ibce(m) == 1) 
a(2*i-1,2*j-1)=a(2*i-1,2* 3-1) -browx (m,2*in-1) ; 
a(2*i-1,2*4)=a(2*i-1,2*j) -browx (m,2*in) ; 
a(2*i,2*j-1)=a(2*i,2* 3-1) -browy (m,2*in-1); 
a(2*i,2*4j)=a(2*i,2* 3) -browy (m,2*in) ; 
btx=btx-arowx (m, 2*in-1) *u(j) -arowx(m,2*in) *v(j); 
bty=bty-arowy (m, 2*in-1) *u(j) -arowy(m,2*in) *v(j); 


end; 


% IF STRESSES ARE PRESCRIBED OVER THE ELEMENT, STORE THE A MATRIX 


r 


% COEFFICIENTS, EXCEPT AT A CORNER NODE WHERE PRESCRIB 


(D DISPLACEMENTS 


% HAVE BEEN IMPOSED. 
if (ibce(m) == 2) 


btx=btx+ (browx (m, 2*in-1) *tmx (m, in) tbrowx (m,2*in) *tmy(m,in) ) / 


bty=bty+ (browy (m, 2*in-1) *tmx (m, in) +browy (m,2*in) *tmy(m,in) ) / 
estore; 

if(ibcn(j) == 2) 
a(2*i-1,2*j-1)=a(2*i-1,2*j-1)+arowx(m,2*in-1); 
a(2*i-1,2*4)=a(2*i-1,2* 3) +arowx(m,2*in); 
a(2*i,2*4j-1)=a(2*i,2*j3-1) +arowy(m,2*in-1); 


a(2*i,2*4j)=a(2*i,2* 3) +arowy (m,2*in) ; 
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end; 
if (ibcen(j) == 1) 


btx=btx-arowx (m, 2*in-1) *u(j) -arowx(m,2*in) *v(j); 


bty=bty-arowy (m, 2*in-1) *u(j) -arowy(m,2*in) *v(j); 
end; 


end; 


end; % each element node in turn 


end; % each element in turn 


% STORE THE B*T COEFFICIENTS AS EXTENSIONS OF MATRIX A. 


a(2*i-1,jmax)=btx; 
a(2*i,jmax)=bty; 


end; % take each node in turn as P 


S$ APPLY NODAL POINT DISPLACEMENT CONSTRAINTS. 
if (nbcpc > 0) 


for ibcpc=1l:nbcpc; % each point constraint in turn 


irow=2* (nodepc (ibcpc) -1) +tidirpe (ibcpc) ; 


for j=l:jmax; % each column in turn 
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a(irow,j)=0.; 

end; 

a(irow, irow)=1.; 

end; 

end; 

Qo 
fe} 

return; 


end function frmtrx 


Boundary Element Methods for Engineers: 
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% each column in turn 


function shape (varargin) 


% 


SUBPROGRAM TO 


EVALUATI 


CI 


TH 


EIR D 


ERIVATIVES, 


AT THE 


AND STORE 


Ap 


% each point constraint in turn 


pendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


VALU 


ES OF THE SHAPE FUNCTIONS AND 


GAUSS POINTS AND NO 


DI] 


global sd sdl sfl egl ngauss sf zg wgl wg 


% 


FS 


STORE 


APPROPRIATE 


% QUADRATURE 


ol? 


INT 


ngauss=8; 


zg (1) =-0.9602898564; 
zg (2)=-0.7966664774; 
zg (3) =-0.5255324099; 
zg (4) =-0.1834346424; 


zg (5)=-zg(4); 
zg (6)=-zg(3); 
zg(7)=-zg(2); 


zg (8)=-zg(1); 
wg (1) =0.1012285362; 
wg (2) =0.2223810344; 
wg (3) =0.3137066458; 
wg (4) =0.3626837833; 
wg (5) =wg (4); 
wg (6) =wg (3); 
wg (7) =wg (2); 
wg (8) =wg (1); 
egl1 (1) =0.013320244; 


COORDINATI 


EGRATION IN XGL AND CGL. 
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en 


FOR LOGARITHMIC FUNCTION 


Boundary Element Methods for Engineers: 
Part 2: Plane Elastic Problems 


egl (3) 
egl (4) 
egl (5) 
egl (6) 
egl (7) 
egl (8) 
wgl (1) 
wgl (2) 
wgl (3) 
wgl (4) 
wgl (5) 
wgl (6) 
wgl (7) 
wgl (8) 


.197871029; 
- 354153994; 
-529458575; 
. 701814530; 
- 849379320; 
- 953326450; 
.164416605; 
-237525610; 
-226841984; 
-175754079; 
112924030; 
-057872211; 
-020979074; 
-003686407; 


% 


% 


NORMAL GAUSSIAN QUADRATURE. 


fo) 


(o) 


for igauss=l:ngauss; 


zeta=zg(igauss) ; 


sf(1,igauss)=0.5*zeta* (zeta-1.); 


sf(2,igauss)=1.-zeta%2; 
uss)=0.5*zeta* (zetat+l.); 


53 


sf(3,iga 


sd(1,igauss)=zeta-0. 


sd(2,igauss)=-2.*ze 


la; 


sd(3,igauss) =zetat+0.5; 


% 


end; each gauss point in turn 


% 


FOUR CASES OF LOGARITHMIC GAUSSIAN Q 


Appendix E: Matlab Version of Quadratic Boundary 
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each gauss point in turn 


= 
a 


r 


UADRATURE TO CONSIDER. 


GI 


ENT FROM FIRST TO THIRD NODE. 


Gl 


NT FROM SECOND TO THIRD NODE. 


Gl 


INT FROM SECOND TO FIRST NODE. 


1 


TO FIRST NODE. 


ENT FROM THIRD 


& IC=1 - INTEGRATION OVER WHOLE ELEM 
& IC=2 - INTEGRATION OVER HALF ELEME 
% IC=3 - INTEGRATION OVER HALF ELEME 
& Ic=4 - INTEGRATION OVER WHOLE ELEM 
for ic=1:4; % each case in turn 


fo) 


(o) 


for igauss=l1:ngauss; 


eta=egl (igauss) ; 


if(ic == 1) 
zeta=2.*eta-1.; 

end; 

if(ic == 2) 


zeta=eta; 
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end; 

if(ic == 3) 
zeta=—eta; 

end; 

if(ic == 4) 
zeta=1.-2.*eta; 

end; 


sfl(ic,1,igauss)=0.5*zeta* (zeta-1.); 


sfl(ic,2,igauss)=1.-zeta%2; 


sfl(ic,3,igauss)=0.5*zeta* (zetatl.); 


sdl(ic,1,igauss)=zeta-0.5; 

sdl(ic,2,igauss)=-2.*zeta; 

sdl(ic,3,igauss)=zetat+0.5; 
end; % each gauss point in turn 


end; % each case in turn 


% 


% SHAPE FUNCTION DERIVATIVES AT THE NODES, STORED AS THOUGH THEY 


% ARE FOR GAUSS POINTS NUMBERED 11, 12 AND 13. 


for igauss=11:13; % each element node in turn 


if(igauss == 11) 


THs @DOOK ts propucep with Wext® 
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zeta=-1.; 
end; 
if(igauss == 12) 
zeta=0.; 
end; 
if(igauss == 13) 
zeta=1.; 
end; 
sd(1,igauss) =zeta-0.5; 
sd(2,igauss)=-2.*zeta; 
sd(3,igauss)=zetat+0.5; 
end; % each element node in turn 
return; 


end function shape 


function intker (i,m) 


% 


% SUBPROGRAM TO INTEGRATE THE KERNEL PRODUCTS FOR A PARTICULAR FORCI 


EA 


E 


% POINT P (INDICATED BY NODE NUMBER I) OVER A PARTICULAR ELEMENT 


% (INDICATED BY ARGUMENT M) BY GAUSSIAN QUADRATURE. 


global node jacob wgl browy browx sfl ngauss bk2yy wg bk2yx bk2xy 
bk2xx sf bkyy bkyx bkxy bkxx akyy arowy akyx akxy arowx akxx 


ungy ungx ynode xnode e nu 


% 


\ 


%& CONSTANT IN DISPLACEMENT KERNELS MULTIPLYING THE LOGARITHMIC TERM. 


cl=(1.+nu)* (3.-nu)/(4.*pix*e); 


% 


% COORDINATES OF POINT P. 


xp=xnode (i); 


yp=ynode (i); 


% 


%& SET UP THE ELEMENT NODE LOOP. 


r 


for in=1:3; % each element node in turn 


j=node(m,in); 


oP 
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& IF P IS NOT THE CURRENT ELEMENT NODE, USE NORMAL GAUSSIAN QUADRATURE. 
if(a = 4) 
for igauss=l:ngauss; % each gauss point in turn 

% EVALUATE JACOBIAN AND UNIT NORMAL VECTOR COMPONENTS, ALSO THE KERNELS 
& AT THE PARTICULAR GAUSS POINT FOR NORMAL QUADRATURE OVER THE WHOLE 
% ELEMENT 

it=1; 

ic=1; 

jacobi(m,igauss,it,ic); 


k 


oe 


foe) 


Ss 
a 
a 


a 
a 
b 
b 
b 


b 


end 
end; 

& IF P 

& IS R 
if (i 

% P AT 
if ( 

S$ TERM 

£ 


ernel (xp, yp,m,igauss,ungx,ungy) ; 


ACCUMULATE THE INTEGRALS. 


fn=sf (in, igauss) ; 

rowx (m,2*in-1) =arowx (m, 2*in-1)+wg (igauss) *akxx*sfn*jacob; 
rowx (m, 2*in) =arowx (m, 2*in) +wg (igauss) *akxy*sfn*jacob; 
rowy (m,2*in-1)=arowy (m, 2*in-1)+wg(igauss) *akyx*sfn*jacob; 
rowy (m, 2*in) =arowy (m, 2*in) +wg (igauss) *akyy*sfn*jacob; 
rowx (m,2*in-1) =browx (m, 2*in-1)+wg (igauss) *bkxx*sfn*jacob; 
rowx (m, 2*in) =browx (m, 2*in) +wg (igauss) *bkxy*sfn*jacob; 
rowy (m,2*in-1) =browy (m, 2*in-1)+wg (igauss) *bkyx*sfn*jacob; 
rowy (m,2*in) =browy (m, 2*in) +wg (igauss) *bkyy*sfn* jacob; 


; % each gauss point in turn 


Gl 


IS THE CURRENT ELEMENT NODE, SOME LOGARITHMIC QUADRATUR 


EQUIRED. 


=) 


THE FIRST NODE OF THE ELEMENT. 


T. 


in ==) 1.) 


S INVOLVING NORMAL QUADRATURE. 

or igauss=l:ngauss; % each gauss point in turn 
1Lt=1; 

ic=1; 


jacobi(m,igauss,it,ic); 
kern2 (m, in, igauss) ; 
sfn=sf(in,igauss); 


browx (m, 2*in-1) =browx (m,2*in-1) +wg (igauss) *bk2xx*sfn*jacob; 
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browx (m, 2* in) =browx (m,2*in) twg (igauss) *bk2xy*sfn*jacob; 
browy (m, 2*in-1) =browy (m, 2*in-1)+wg (igauss) *bk2yx*sfn* jacob; 
browy (m, 2* in) =browy (m, 2*in) +wg (igauss) *bk2yy*sfn*jacob; 


end; % each gauss point in turn 


oe 


ERMS INVOLVING LOGARITHMIC QUADRATURE. 


oe 
=| 


for igauss=l:ngauss; % each gauss point in turn 
it=2; 
ic=1; 
jacobi(m,igauss,it,ic); 
sfn=sfl(ic,in,igauss) ; 
dzde=2.; 
browx (m, 2*in-1) =browx (m,2*in-1) +wgl (igauss) *cl*sfn*jacob*dzde; 
browy (m, 2* in) =browy (m, 2*in) +wgl (igauss) *cl*sfn*jacob*dzde; 


end; % each gauss point in turn 


end; 
% P AT THE SECOND NODE OF THE ELEMENT. 
if(in == 2) 
“S\ 
¥ 
\ \\ 
| 
| 
VA 
Deloi 
eloitte. 
Discover the truth at www.deloitte.ca/careers © Deloitte & Touche LLP and affiliated entities. 
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foe) 
= 


ERMS INVOLVING NORMAL QUADRATURE. 


fo) 


for igauss=l:ngauss; % each gauss point in turn 
LHL; 
ic=1; 
jacobi(m,igauss,it,ic); 
kern2 (m, in, igauss) ; 
sfn=sf(in,igauss); 
browx (m, 2*in-1) =browx (m,2*in-1) +wg (igauss) *bk2xx*sfn*jacob; 
browx (m, 2* in) =browx (m, 2*in) twg (igauss) *bk2xy*sfn*jacob; 


browy (m, 2*in-1) =browy (m,2*in-1) +wg (igauss) *bk2yx*sfn*jacob; 


browy (m, 2* in) =browy (m, 2*in) twg (igauss) *bk2yy*sfn*jacob; 


fo) 


end; % each gauss point in turn 


% TERMS INVOLVING LOGARITHMIC QUADRATURE. 


fo) 


for igauss=l:ngauss; % each gauss point in turn 
LeE=27 
ic=2; 
jacobi(m,igauss,it,ic); 
sfn=sfl(ic,in,igauss) ; 
dzde=1.; 
browx (m, 2*in-1) =browx (m,2*in-1) +wgl (igauss) *cl*sfn*jacob*dzde; 
browy (m, 2* in) =browy (m, 2*in) +wgl (igauss) *cl*sfn*jacob*dzde; 
ic=3; 
jacobi(m,igauss,it,ic); 
sfn=sfl(ic,in,igauss) ; 
dzde=1.; 
browx (m, 2*in-1) =browx (m,2*in-1) +wgl (igauss) *cl*sfn*jacob*dzde; 
browy (m, 2* in) =browy (m, 2*in) twgl (igauss) *cl*sfn*jacob*dzde; 
end; % each gauss point in turn 


end; 


ol? 


r 


ol? 


P AT THE THIRD NODE OF THE ELEMENT. 


if(in == 3) 


oP 


Gl 


foe) 
= 


ERMS INVOLVING NORMAL QUADRATURE. 


fo) 


for igauss=l:ngauss; % each gauss point in turn 
LeE=1F 
ic=1; 
jacobi(m,igauss,it,ic); 
kern2 (m, in, igauss) ; 
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sfn=sf(in,igauss); 

browx (m, 2*in-1) =browx (m,2*in-1) +wg (igauss) *bk2xx*sfn*jacob; 
browx (m, 2* in) =browx (m,2*in) +wg (igauss) *bk2xy*sfn*jacob; 
browy (m, 2*in-1) =browy (m, 2*in-1) +wg (igauss) *bk2yx*sfn*jacob; 
browy (m, 2* in) =browy (m, 2*in) +wg (igauss) *bk2yy*sfn*jacob; 


fo) 


end; % each gauss point in turn 


% TERMS INVOLVING LOGARITHMIC QUADRATURE. 
for igauss=1l:ngauss; % each gauss point in turn 
it=2; 
ic=4; 
jacobi(m,igauss,it,ic); 
sfn=sfl(ic,in,igauss) ; 
dzde=2.; 
browx (m, 2*in-1) =browx (m,2*in-1)+wgl (igauss) *cl*sfn*jacob*dzde; 
browy (m, 2* in) =browy (m, 2*in) twgl (igauss) *cl*sfn*jacob*dzde; 
end; % each gauss point in turn; 
end; 
end; 
end; % each element node in turn 


return; 


end Sfunction intker 


function jacobi(m,igauss,it,ic) 


% 


fo) 


SUBPROGRAM TO EVALUATE THE JACOBIAN AND THE COMPONENTS OF THE UNIT 


r 


= 


% NORMAL VECTOR AT A PARTICULAR GAUSS POINT. 


r 


% M INDICATES THE ELEMENT NUMBER. 


% IGAUSS INDICATES THE 


zr 


GAUSS POINT NUMB 


ea 
w 


% IT INDICATES THE TYPE OF QUADRATURE. 


% IT=1 - NORMAL GAUSSIAN QUADRATURE. 
% IT=2 - LOGARITHMIC GAUSSIAN QUADRATURE. 


% IC INDICATES THE CASE NUMBER FOR LOGARITHMIC GAUSSIAN QUADRATUR 


ea 


% AS DEFINED IN SUBROUTINE SHAPE. 


oP 


global jacob ungy ungx node ynode xnode sdl sd 
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% CALCULATE THE DERIVATIVES OF THE 
% THE LOCAL INTRINSIC COORDINATE. 
dx=0.; 

dy=0.; 


for in=1:3; 
if(it == 1) 
sfnd=sd(in,igauss) ; 
end; 
if(it == 2) 
sfnd=sdl (ic,in,igauss) ; 
end; 
j=node(m,in); 
dx=dx+sfnd*xnode (j); 
dy=dy+sfnd*ynode (4); 
end; % each element node in turn 


% 


% COMPONENTS OF THE 


ungx=dy; 
ungy=-dx; 
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GLOBAL COORDINATES WITH RESPECT TO 


% each element node in turn 


GAUSS POINT. 


> 
accenture 


High performance. Delivered. 


173 


© 
Ny 
= 
w 
> 
2 
Q 
oO 
= 
= 
= 
8 


Boundary Element Methods for Engineers: 
Part 2: Plane Elastic Problems 


% 


% 


JACOBIAN OF THE 


r 


COORDINATE 


r 


jacob=sqrt (ungx*2+ungy*2) ; 


% 


Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


TRANSFORMATION. 


%& SCALE THE VECTOR COMPONENTS TO GIVE THE UNIT NORMAL VECTOR. 
ungx=ungx/jacob; 

ungy=ungy/jacob; 

return; 

end function jacobi 

function kernel (xp, yp,m, igauss, unx, uny) 

% SUBPROGRAM TO EVALUATE THE KERNELS WHEN P IS NOT THE CURRENT 
% ELEMENT NODE. 

% XP, YP INDICATE THE GLOBAL COORDINATES OF POINT P. 

% M INDICATES THE ELEMENT NUMBER. 

% IGAUSS INDICATES THE NUMBER OF THE GAUSS POINT, Q. 

global bkyy bkxy bkyx bkxx akyy akyx akxy akxx nu node ynode xnode sf 


% 


% 


COORDINATES OF GAUSS POINT QO 


xq=0.; 
yq=0.; 

for in=1:3; % each element node 
sfn=sf(in,igauss); 
j=node(m,in); 


xg=xqtsfn*xnode (4); 


yq=yqtsfn*ynode (j); 


% 


end; 


% 


in turn 


each element node in turn 


% 


COMPONENTS AND MAGNITUDE 


I 


UX=XQ-Xp; 
LY=Yqd-YP; 
rsq=erx*2+ry*2; 
r=sqrt(rsq); 
LX=Cx/T; 


ry=ry/v; 
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% 


NORMAL TO TH 


Appendix E: Matlab Version of Quadratic Boundary 
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= 
py 


zr 


BOUNDARY AT Q. 


%& RATE OF CHANGE OF R WITH THE 
drdn=rx*unxtry*uny; 
% PARAMETERS IN THE KERNEL FUNCTIONS. 


cl=-1./(4.*pi*r); 
c2=1.-nu; 
c3=2.*(1.+nu); 
c4=(1.+nu)%*2/(4.*pix*e); 


Go 


(3.-nu)*log(1./r)/(1.+nu) ; 


% 


z E 


EVALUATE 


G 


THE 


K 


ERN 


ELS. 


akxx=cl* (c2+c3*rx*rx) *drdn; 
terml=c3*rx*ry*drdn; 
term2=c2* (rx*uny-ry*unx) ; 


kxy=cl* (terml-term2) ; 


kyx=cl* (terml+term2) ; 
kyy=cl* (c2+c3*ry*ry) *drdn; 
kxx=C4* (c5+rx* rx) ; 
kxy=Cc4*rx* ry; 


kyx=bkxy; 


kyy=c4* (cotry*ry) ; 


foe) 


return; 


end %$function kernel 


function kern2(m,in,igauss) 

% SUBPROGRAM TO EVALUATE THE NON-SINGULAR LOGARITHMIC TERM IN THE 
% SECOND KERNEL WHEN P IS THE CURRENT ELEMENT NODE. 

% M INDICATES THE ELEMENT NUMBER. 

% IN INDICATES THE NUMBER OF THE ELEMENT NODE FORMING P. 

% IGAUSS INDICATES THE GAUSS POINT NUMBER. 

global bk2yy bk2xy bk2yx bk2xx nu ynode xnode zg node sf 


% 


% 


COORDINAT! 


ES OF GAUSS POINT 0. 
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xq=0.; 

yq=0.; 

for ic=1:3; % each element node in turn 
sfn=sf(ic,igauss); 
j=node(m,ic); 
xg=xqtsfn*xnode (4) ; 
yq=yqtsfn*ynode (4); 


fo) 


end; % each element node in turn 


% 


% ELEMENT NODE NUMBERS. 


il=node(m,1); 
i2=node (m, 2); 
i3=node (m, 3); 


% 


% EVALUATE THE INTRINSIC COORDINATE. 


zeta=zg(igauss) ; 


ol? 


= 


6 P AT THE FIRST NODE OF THE ELEMENT. 
if(in == 1) 


xcomp=(zeta-2.) *xnode(il)+2.* (1.-zeta) *xnode (i2) +zeta*xnode (i3) ; 


The Wake 
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ycomp=(zeta-2.) *ynode(il)+2.* (1.- 


rfin=sqrt (xcomp*2+ycomp%2) ; 


Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


zeta) *ynode (i2)+zeta*ynode (13) ; 


ECOND NODI 


OF TH 


EM 


ENT. 


xcomp=-0.5* (zeta-1.) *xnode(il)+ze 
ycomp=-0.5* (zeta-1.) *ynode(il)+ze 


rfin=sqrt (xcomp*2+ycomp%2) ; 


iy 
Bi 


= 
py 


P AT TH! 


THIRD NOD! 


r 


OF THE 


EM 


) 


xcomp=-zeta*xnode (il)+2.* (1.+zeta 


if(in == 


ycomp=-zeta*ynode (i1l)+2.*(1.+zeta 
rfn=sqrt (xcomp*2+ycomp%2) ; 
i1=13; 


end; 


6 


COMPON 


ENTS AND MAGNITUDI 


I 


OF THI 


G 


py py 


xp=xnode (i); 
yp=ynode (i); 
UX=XQ-Xp; 
LY=Yqd-YPr 
rsq=rx*2t+ry%*2; 
r=sqrt(rsq); 
LX=Cx/X; 


ry=ry/vr; 


ol? 


ET EL 


G 


ERS IN THE 
u) *2/ (4. *pi*e); 


K 


ERN 


FUNCTIO 


u) *log(1./rfn)/(1.+nu); 


G 


EVALUATE 


G 


THE 


K 


ERN 


oP 


ELS. 


bk2xx=c4* (c5t+rx*rx) ; 

bk2xy=c4*rx* ry; 

bk2yx=bk2xy; 

bk2yy=c4* (co+ry*ry); 
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ta*xnode (i2)-0.5* (zetat+1l.) *xnode(i3); 


ta* ynode (i2)-0.5* (zeta+1.) *ynode (i3) ; 


ENT. 


) *xnode (i2)-(zetat2.) *xnode (i3); 
) *ynode (i2)-(zetat+2.) *ynode (i3); 
RADIUS VECTOR FROM P TO Q. 


NS. 


77 


Boundary Element Methods for Engineers: Appendix E: Matlab Version of Quadratic Boundary 
Part 2: Plane Elastic Problems Element Program for Plane Elastic Problems 


return; 


end function kern2 


function output (varargin) 


% 


= 


% SUBPROGRAM TO WRITE OUT THE NODAL POINT VALUES OF DISPLACEMENTS 


% AND ELEMENT STRESSES AND COMPUTE THE FORCES ON THE BOUNDARY SEGMENTS. 


% 


global nsegtot fyseg fxseg node maxl tmy jacob wg tmx sf ngauss 


isegelem nel sige sigsn sigss signn nu e unmx unmy v sd u ty tx 


ibce estore nnp uv ibcn fid6 


% ARRANGE FOR U, V AND TX, TY TO CONTAIN THE NODAL DISPLACEMENTS 


% AND TRACTIONS, RESPECTIVELY. 


for i=l:nnp; % each node in turn 
if(ibcn(i) == 1) 
tx (Li) =uv(2*i-1); 
ty (1) =uv(2*i); 


end; 


u(i)=uv(2*i-1); 
v(i)=uv(2*i); 
end; 


end; % each node in turn 


% 


% HEADING FOR OUTPUT OF NODAL DISPLACEMENTS AND TRACTIONS. 


fprintf (fid6, ['\n', 'NODAL POINT DISPLACEMENTS AND TRACTIONS','\n','\n', 


: I U Vv TX ry, "An! ]) 3 
e=estore; 


for i=l:nnp; % each node in turn 


% 


% REMOVE THE SCALING. 


u(i)=u(i) *maxl; 
v(i)=v(i)*maxl; 
tx (1) =tx (i) *estore; 
ty (i) =ty (i) *estore; 
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% OUTPUT THE NODAL VALUES OF DISPLACEMENTS AND TRACTIONS. 


fprintf (fid6, ['S4i','S15.6e', 'S15.6e', 'S15.6e', 'S15.6e', 


"\n'j,i,u(i),v(i),tx(i),ty(i)); 


end; % each node in turn 


% 


% HEADING FOR OUTPUT OF ELEMENT STRESSES. 


fprintf (fide, ['\n', 'STRESSES AT THE NODES OF THE ELEMENTS’, 


Nat pt Nn M IN SIGNN SIGSS SIGSN', 


! SIGE', '\n']); 


so 
() 


% NORMAL AND SHEAR STRESSES AT THE NODES OF THE ELEMENTS. 


for m=l:nel; % each element in turn 
for in=1:3; % each element node in turn 

if (ibce(m) ~= 2) 
j=node(m,in); 
tmx (m, in) =tx (4); 
tmy (m, in)=ty(j); 
signn (m, in) =tmx (m, in) *unmx (m,in)+tmy(m,in) *unmy(m,in); 
sigsn(m, in) =-tmx (m, in) *unmy (m,in)+tmy(m,in) *unmx(m, in) ; 


end; 
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% 


% DIRECT STRESS ALONG THE BOUNDARY SURFACE. 


r 


EA 


icl=node(m,1); 

ic2=node (m, 2); 

ic3=node(m,3); 

igauss=in+10; 

dudz=sd(1,igauss) *u(icl)+sd(2,igauss) *u(ic2)+sd(3,igauss) *u(ic3) ; 
dvdz=sd(1,igauss) *v(icl)+sd(2,igauss) *v(ic2)+sd(3,igauss) *v(ic3) ; 
LEH LF 

ic=1; 

jacobi(m,10+in,it,ic); 
ess=(-dudz*unmy (m, in) tdvdz*unmx (m, in) ) / (jacob*max1) ; 


sigss(m,in)=e*esstnu*signn (m, in); 


ol? 


oP 


VON MISES EQUIVALENT STRESS. 


sige (m,in)=sqrt(signn(m, in) *2+sigss(m,in) *2 


-signn(m,in)*sigss(m,in)+3.*sigsn(m,in)%2) ; 


foe) 


oP 


OUTPUT THE STRESSES AT THE NODES OF THE ELEMENTS. 
fprintf (fid6é, ['S4i','S4i','%S15.6e', 'S15.6e', '$15.6e','S15.6e', 


'\n'j,m,in,signn(m,in),sigss(m,in),sigsn(m,in),sige(m,in)); 
end; % each element node in turn 
fe) 


end; % each element in turn 


% 


% COMPUTE THE FORCES ON TH 


GI 


BOUNDARY SEGMENTS. 


fo) 


for iseg=l:nsegtot; % each segment in turn 
fxseg(iseg)=0.; 
fyseg(iseg)=0.; 

end; % each segment in turn 


for m=l:nel; % each element in turn 


iseg=isegelem(m) ; 


oe 


% APPLY GAUSSIAN QUADRATURE. 


Gl 


fxelem=0.; 
fyelem=0.; 
for in=1:3; % each element node in turn 
for igauss=l:ngauss; % each gauss point in turn 
sfn=sf(in,igauss); 
it=1; 
ic=1; 
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jacobi(m,igauss,it,ic); 

fxelem=fxelemt+wg (igauss) *sfn*jacob*tmx (m, in) *maxl; 

fyelem=fyelemtwg (igauss) *sfn*jacob*tmy(m, in) *maxl; 
end; % each gauss point in turn 


end; % each element node in turn 


% 


\ 


6 ACCUMULAT 


= 


THE FORCES ON THE SEGMENT. 


fxseg (iseg) =fxseg(iseg)+fxelem; 


fyseg (iseg) =fyseg(iseg)+fyelem; 


end; % each element in turn 


% 


% OUTPUT TH 


zr 


SEGMENT FORCE RESULTS. 


r 


fprintf (fid6é, ['\n','FORCES ACTING ON THE BOUNDARY SEGMENTS', 


"\n', '\n', 'SEGMENT FX FY', '\n']); 


for iseg=1l:nsegtot; 
fprintf (fid6é, ['S5i','%14.5e','S14.5e','\n'J, 
iseg,fxseg(iseg),fyseg(iseg) ); 
end; 
return; 


end function output 


function [x,iflag]=elimin(a,neqn) 


% 


oe 
x 


SUBROUTINE FOR SOLVING SIMULTANEOUS LINEAR EQUATIONS BY GAUSSIAN 


% ELIMINATION WITH PARTIAL PIVOTING. 


foe) 


INITIALIZ! 


Gl 


ILL-CONDITIONING FLAG. 
iflag=0; 
for i=l:neqn; x(i)=0.; end; 


% 


% SCALE EACH EQUATION TO HAVE A MAXIMUM COEFFICIENT MAGNITUDE OF UNITY. 


jJmax=negqnt1; 
for i=l:neqn; % each equation in turn 
amax=0.; 
ie) 


for j=l:neqn; % search for maximum 


absa=abs (a(i,j)); 


Download free eBooks at bookboon.com 


181 


Boundary Element Methods for Engineers: Appendix E: Matlab Version of Quadratic Boundary 
Part 2: Plane Elastic Problems Element Program for Plane Elastic Problems 


if(absa > amax) 
amax=absa; 
end; 


end; % search for maximum 


oe 


fo) 


for j=l:jmax; % scale coefficients 
a(i,j)=a(i,j) /amax; 
end; % scale coefficients 


end; % each equation in turn 


% 


% COMMENCE ELIMINATION PROCESS. 


fo) 


for k=l:neqn-1; % eliminate each variable in turn 


oP 
1p) 


EARCH LEADING COLUMN OF THE COEFFICIENT MATRIX FROM THE DIAGONAL 


% DOWNWARDS FOR THE LARGEST VALUE. 


= 


imax=k; 
for i=k+l:neqn; % search for largest value 
if(abs(a(i,k)) > abs(a(imax,k))) 


imax=i; 
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end; 


end; % search for largest value 


oe 


THE LARGEST COEFFICIENT 


r 


r 


fo) 


IF NECESSARY, INTERCHANG 


EQUATIONS TO MAK 


BECOME THE PIVOTAL COEFFICIENT. 


oe 


if (imax ~= k) 
for j=k:jmax; % interchange coefficients 
atemp=a(k,j); 
a(k,j)=a(imax, J); 
a (imax, j)=atemp; 
end; % interchange coefficients 


end; 


= 


% ELIMINATE X(K) FROM EQUATIONS (K+1) TO NEON, FIRST TESTING FOR 


%& EXCESSIVELY SMALL PIVOTAL COEFFICIENT (ASSOCIATED WITH A SINGULAR 


% OR VERY ILL-CONDITIONED MATRIX). 


if (abs(a(k,k)) < 1.0e-5) 
iflag=1; 
return; 
end; 
for i=k+l:neqn; % each of remaining equations 
fact=a(i,k)/a(k,k); 
for j=k:jmax; % modify coefficients 
a(i,j)=a(i,j)-fact*a(k,j); 
end; % modify coefficient 
end; % each of remaining equations 


% 


end; % eliminate each variable in turn 


% 


% SOLVE THE EQUATIONS BY BACK SUBSTITUTION, FIRST TESTING 


% FOR AN EXCESSIVELY SMALL LAST DIAGONAL COEFFICIENT. 


if (abs (a(neqn,negn)) < 1.0e-5) 
iflag=1; 
return; 

end; 

x (neqn) =a (neqn, jmax) /a(neqn,neqn) ; 

for i=negqn-1:-1:1; % then each unknown in turn backwards 
sum=a (i, j}max) ; 


for j=i+l:neqn; % sum products 


sum=sum-a (i,j) *x(j); 
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fo) 


end; % sum products 
x(i)=sum/a(i,i); 

end; 

return; 


end function elimin 


Appendix E: Matlab Version of Quadratic Boundary 
Element Program for Plane Elastic Problems 


% then each unknown in turn backwards 


The equivalent of the Fortran SHAREDDDATA2EQ is shareddata2eq, stored as the file shareddata2eq.m, 


is as follows. 


[oume) 


%%6%S--- module shareddata2Zeq; 


% module STORING SHARED DATA. 


global xeend; if isempty(xeend), xeend=zeros(1,260); end; 
global yeend; if isempty(yeend), yeend=zeros(1,260); end; 
global xnode; if isempty(xnode), xnode=zeros(1,500); end; 
global ynode; if isempty(ynode), ynode=zeros(1,500); end; 
global xsend; if isempty(xsend), xsend=zeros(1,250); end; 
global ysend; if isempty(ysend), ysend=zeros(1,250); end; 
global unmx; if isempty(unmx), unmx=zeros(250,3); end; 
global unmy; if isempty(unmy), unmy=zeros(250,3); end; 
global angstore; if isempty(angstore), angstore=zeros(1,260); end; 
global maxl; if isempty(maxl), maxl=0; end; 

global a; if isempty(a), a=zeros(1000,1001); end; 

global uv; if isempty(uv), uv=zeros(1,1000); end; 

global u; if isempty(u), u=zeros(1,500); end; 

global v; if isempty(v), v=zeros(1,500); end; 

global useg; if isempty(useg), useg=zeros(1,250); end; 
global vseg; if isempty(vseg), vseg=zeros(1,250); end; 
global arowx; if isempty(arowx), arowx=zeros(250,6); end; 
global arowy; if isempty(arowy), arowy=zeros(250,6); end; 
global browx; if isempty(browx), browx=zeros(250,6); end; 
global browy; if isempty(browy), browy=zeros(250,6); end; 
global zg; if isempty(zg), zg=zeros(1,8); end; 

global wg; if isempty(wg), wg=zeros(1,8); end; 

global egl; if isempty(egl), egl=zeros(1,8); end; 

global wgl; if isempty(wgl), wgl=zeros(1,8); end; 

global jacob; if isempty(jacob), jacob=0; end; 

global ungx; if isempty(ungx), ungx=0; end; 
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global ungy; if isempty(ungy), ungy=0; end; 

global signnseg; if isempty(signnseg), signnseg=zeros(1,250); end; 
global sigsnseg; if isempty(sigsnseg), sigsnseg=zeros(1,250); end; 
global uelem; if isempty(uelem), uelem=zeros(1,250); end; 

global velem; if isempty(velem), velem=zeros(1,250); end; 

global signn; if isempty(signn), signn=zeros(250,3); end; 

global sigss; if isempty(sigss), sigss=zeros(250,3); end; 

global sigsn; if isempty(sigsn), sSigsn=zeros(250,3); end; 

global sige; if isempty(sige), sige=zeros(250,3); end; 

global tx; if isempty(tx), tx=zeros(1,500); end; 

global ty; if isempty(ty), ty=zeros(1,500); end; 

global tmx; if isempty(tmx), tmx=zeros(250,3); end; 

global tmy; if isempty(tmy), tmy=zeros(250,3); end; 

global fxseg; if isempty(fxseg), fxseg=zeros(1,250); end; 

global fyseg; if isempty(fyseg), fyseg=zeros(1,250); end; 

global sf; if isempty(sf), sf=zeros(3,8); end; 

global sd; if isempty(sd), sd=zeros(3,13); end; 

global sfl; if isempty(sfl), sfl=zeros(4,3,8); end; 

global sdl; if isempty(sdl), sdl=zeros(4,3,8); end; 

global e; if isempty(e), e=0; end; 
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global nu; if isempty(nu), nu=0; end; 

global estore; if isempty(estore), estore=0; end; 

global akxx; if isempty(akxx), akxx=0; end; 

global akxy; if isempty(akxy), akxy=0; end; 

global akyx; if isempty(akyx), akyx=0; end; 

global akyy; if isempty(akyy), akyy=0; end; 

global bkxx; if isempty(bkxx), bkxx=0; end; 

global bkxy; if isempty(bkxy), bkxy=0; end; 

global bkyx; if isempty(bkyx), bkyx=0; end; 

global bkyy; if isempty(bkyy), bkyy=0; end; 

global bk2xx; if isempty(bk2xx), bk2xx=0; end; 

global bk2xy; if isempty(bk2xy), bk2xy=0; end; 

global bk2yx; if isempty(bk2yx), bk2yx=0; end; 

global bk2yy; if isempty(bk2yy), bk2yy=0; end; 

global nel; if isempty(nel), nel=0; end; 

global nnp; if isempty(nnp), nnp=0; end; 

global maxnel; if isempty(maxnel), maxnel=0; end; 

global maxnnp; if isempty(maxnnp), maxnnp=0; end; 

global maxnb; if isempty(maxnb), maxnb=0; end; 

global neend; if isempty(neend), neend=0; end; 

global ngauss; if isempty(ngauss), ngauss=0; end; 

global node; if isempty(node), node=zeros(250,3); end; 
global ml; if isempty(ml1), ml=zeros(1,250); end; 

global m3; if isempty(m3), m3=zeros(1,250); end; 

global nbound; if isempty(nbound), nbound=0; end; 

global nsegtot; if isempty(nsegtot), nsegtot=0; end; 
global nelb; if isempty(nelb), nelb=zeros(1,10); end; 
global nsegb; if isempty(nsegb), nsegb=zeros(1,10); end; 
global nbcu; if isempty(nbcu), nbcu=0; end; 

global nbcs; if isempty(nbcs), nbcs=0; end; 

global nbcm; if isempty(nbcm), nbcm=0; end; 

global nbct; if isempty(nbct), nbct=0; end; 

global ibce; if isempty(ibce), ibce=zeros(1,250); end; 
global ibcn; if isempty(ibcn), ibcn=zeros(1,500); end; 
global isegbc; if isempty(isegbc), isegbc=zeros(1,250); end; 
global isegend; if isempty(isegend), isegend=zeros(1,250); end; 
global isegelem; if isempty(isegelem), isegelem=zeros(1,250); end; 
global mfirst; if isempty(mfirst), mfirst=zeros(1,250); end; 
global mlast; if isempty(mlast), mlast=zeros(1,250); end; 
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Part 2: Plane Elastic Problems Element Program for Plane Elastic Problems 
nbcpc; isempty(nbcpc), nbcpc=0; ; 
maxnpc; isempty(maxnpc), maxnpc=0; 1; 
nodepc; isempty(nodepc), nodepc=zeros (1,20); ; 
idirpc; isempty(idirpc), idirpc=zeros(1,20); end; 
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Solutions to Problems — Part Il 


Numerical answers are given, followed where appropriate by detailed notes on the solution procedure 


and comments on the results obtained. 


Chapter 5 


5.1 Indistinguishable when v = 0 (v* = v and E* = E). Oflittle practical significance because no 


known materials have Poisson's ratios of zero, although some have values close to zero. 


52 An incompressible material has a Poisson's ratio of v = 1/2. In plane strain v* = 1 and 


E* = 4E/3. The typical displacement and traction kernel functions become 


Use (Bq) = [In (=) + AF, | 


Be AK 
Uyy (p,q) = tne *'Y 


1,., dr 
Tyx (P,Q) = Se a 
1, , dr 
Try (p, 9) = ae "y in 


Unlike displacement-based finite element methods, incompressible plane strain presents no 


difficulties to the present boundary element formulation. 
5.3 tyAy + tyAy = —k, (un, + vi) 
igny — ih, = —k, (on, — un, 
5.4 The body forces per unit volume in Equations 1.6 and 1.7 are X = 0 and Y = —pg, and the 
expressions for stresses satisfy these equations. Strains e,, and e,, are zero, from both the strain 


definitions of Equations 1.2 and 1.3, and in terms of stresses from Equations 1.20 and 1.23. From 


Equations 1.2 


and from Equation 1.21 


1 ‘ pay j 
Cyy = Fr (%y — V*oxx ) ==, (A-yv 2 


The expressions for ey, are identical, and all the relevant equations are satisfied. 
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5.5 


The body forces per unit volume in Equations 1.6 and 1.7 are X = pw*x and Y = 0, and the 
expressions for stresses satisfy these equations. Strains e,, and e,, are zero, from both the strain 
definitions of Equations 1.2 and 1.3, and in terms of stresses from Equations 1.21 and 1.23. From 
Equations 1.2 

pw*x? 


Cx, = OB* (1 — v*?) 


and from Equation 1.20 


pw*x? 


(4 — y* 
apr (EV) 


1 
Cxx = BE (Sy = V* Oxy ) = 


The expressions for e,, are identical, and all the relevant equations are satisfied. 


Chapter 6 


6.1 


6.2 


Answer: maximum hoop stress concentration factor 1.005. The analytical solution for a hole 
in an infinite plate (K — oo) is 1 (Equations 6.8). For a plate width 20 times the hole diameter 
(K = 20), which ignores the square corners of the plate, it is 1.005, identical to four significant 
figures to that computed. But, with compressive radial stress (the pressure) at the surface of the 
hole, the hoop stress is not the largest stress to be considered there. The maximum von Mises 
equivalent stress concentration factor (relative to the applied pressure) is 1.74. Results were 


obtained using 4 elements per side of the outer boundary, 6 elements per quadrant of the hole. 


Answers: at the inner surface, hoop stress 0.231, radial displacement 0.00240. At the outer 
surface, radial stress -0.539, hoop stress -0.231. Stress values are averaged over element end 
nodes and midside nodes: for example, typical hoop stress values at the inner surface are 0.2313 
at end nodes, 0.2309 at midside nodes (a difference which is not significantly reduced by mesh 


refinement). Results were obtained using 4 elements per quadrant on both boundaries. 


The general analytical solutions for stresses in a thick-walled cylinder are given by Equations 
6.5. In this problem the boundary conditions are o,,=—p at r=1, and egg =0 at 


Y = 1%, from which 


A= — ee rt B= Partin? 
(1—-v*)r£+(1+v*)r? ’ (-v)rZ+(14v*)r? 


At the inner surface,r = r, and the hoop stress is 


_ p[(G-v*)K?-(14v*)] 
700 = y+ 4) 
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where, K = 2 = 2,p = 1 and y* = > = 0.42857, giving ogg = 0.2308. 


Ty —vV 


At the outer surface, r = ry and the radial and hoop stresses are 


| ee === 
Orr = Gy REA) = 70-5384, 999 = "Oy = — 0.2308 


Radial displacement at the inner surface can be found from the hoop strain with the aid of 


Equation 6.9 


% 
uy = a oe — V" Ory) 


“= 1098.9, giving u, = 0.002402. 


Where 099 = 0.2308, o, = —1, and E™ = (1-v2) ~ 


6.3 Answer: 3.12. Result obtained using 4 elements per edge of the plate on edges without the notch, 
4 elements per quadrant in the notch, and 10 elements on each of the segments adjacent to the 
notch. No significant change in the result was found when element sizes on these two segments 


were graded to give smaller elements near the notch. 


6.4 Answer: 7.81. Result obtained using 4 elements per outer edge of the plate, 20 per straight side 


of the square hole, and 8 per quadrant of the rounded corners. 
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6.5 


6.6 


6.7 


6.8 


Answer: 0.681, a friction coefficient of more than 0.5 would be required to prevent slipping. 
The friction requirement is investigated by finding the ratio between shear and normal stresses 
(SIGSN and SIGNN) at nodes on the base of the pedestal, giving the minimum coefficient to 
prevent slipping. In this case the value of 0.5 is exceeded at points more than about 0.07 m 
from the central line of symmetry (less than about 0.08 m from the edges of the pedestal base). 
Results were obtained using 4 elements on each of the top and sloping sides of the pedestal, 
and 8 on the base. The boundary conditions are zero displacements on the base, and a normal 
stress of -1 on the top of the pedestal. Such contact problems, involving slippage with friction 
and possible local loss of contact, can be solved by boundary element methods, but are non- 
linear in the sense that the boundary conditions change with the loading, and require more 


sophisticated programming. 


Answers: stiffness 142 MN/m, maximum stress (SIGSS) 190 MN/m/ at the fillets, stress (SIGSS) 
variation over the central 30 mm of length is from 89.3 to 89.7 MN/m’, compared with the 
nominal stress of 90 MN/m7”. A state of plane stress may be assumed. Results were obtained using 
4 elements per straight edge at the ends of the test piece, 4 elements on each of the fillets, and 10 
elements along each side of the central section (to provide more detail there, and to have element 
nodes at the ends of the central 30 mm). Point constraints were applied to the corners at the left 
hand end of the test piece. The stiffness is estimated from the computed relative displacement 
between the nodes at the centres of the two ends: 0.380 x 10~* m. The tensile force in the test 
piece is 60 x 10° x 0.03 x 0.003 = 5400 N (stress times width times thickness), giving a 
stiffness of 5400/0.380 x 10~* N/m or 142 MN/m. 


To apply pure shear, balancing shear stresses must be applied to all four edges of the plate: 
for example,o,, = +1 on BC and DA, 0,, = —1 on CD and AB (Figure 6.7). The maximum 
computed hoop stresses at the hole are 4.03, occurring at four angular positions at 45° to the 
x and y axes. The analytical solution for a hole at the centre of an infinite plate in pure shear 
gives a value of 4 times the applied shear stress. Results were obtained using 4 elements per side 


of the outer boundary, 6 elements per quadrant of the hole. 


Answer: 5.03 (hoop stress concentration factor at both sides of the hole). Rather than create 
the mesh data for an elliptical hole boundary, the same geometric dataas in Problems 6.1 and6.7 can 
used, but with the y coordinates of the corners of the outer boundary doubled in magnitude. Then in 
subprogram MESHQall YNODE values can be scaled by a factor of 0.5 immediately before finding 
the maximum dimension of the domain. This creates the elliptical shape ofhole and brings the outer 
boundary shape back to square. Results were obtained using 4 elements per sideoftheouterboundary, 
6 elements per quadrant of the hole. The exact solution for an elliptical hole with a ratio of major to 


minor axes of 2, under tension in an infinite plate, is 5. 
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6.9 


6.10 


Answers: (a) 3.11 at the outer sides of the holes (y = 0,x = +0.075), 3.06 at the inner sides 
(y = 0,x = +0.025). (b) 2.74 (at x = +0.05,y = +0.025). Results were obtained using 4 
elements per side of the outer boundary, 6 elements per quadrant of each of the holes. An 


example of a solution domain with three boundaries. 


Answer: 745 MN/m’ at the bottom of the groove. The position of the maximum equivalent stress 
is as expected. With a radius ratio of K = 2.5, the hoop stresses at the inner and outer surface 


of the cylinder in the absence of the groove are 


2 2 
O99 =P oD = 166MNi/m MN/m? and %@g = Ga = 45-7 MN/m’ 


These figures are useful for checking the computed results: hoop stresses at the inner surface 
remote from the groove should be reasonably close to these (within about 1-2%, say). The 
choice of mesh in this problem requires some care. Experience shows that 4 elements per 
quadrant should be sufficient on the outer boundary, and perhaps 6 elements per quadrant 
in the semi-circular groove. On the inner cylindrical surface it is important to have elements 
near the groove which are similar in length to the very small elements in the groove. Using 
elements of uniform size to satisfy this requirement leads to an unnecessarily large number of 
elements. A reasonable compromise was found to be 16 elements on each of two nearly semi- 
circular segments (that is, about 8 per quadrant), with an element length ratio (as explained in 
Sections 3.1.3 and 4.1.3) of. In other words, the constant ratio between the lengths of successive 
elements moving away from the groove was 1.3. The total number of elements used was 60. In 
addition to checking inner surface hoop stresses remote from the groove against the analytical 
values above, another useful test is for the tangential stresses (SIGSS) in the elements adjacent 
to the corner between the groove and the inside surface of the cylinder: at the common node 
they should both be close to the applied pressure. The computed hoop stress at the bottom of 
the groove is 678 MN/m’, 4.08 times the hoop stress at the inner surface in the absence of the 
groove. Combined with a radial stress of -120 MN/m’, this gives a von Mises equivalent stress 
of 745 MN/m/?, some 6.2 times the applied pressure in magnitude, also some 3.00 times the 


maximum equivalent stress in the absence of the groove. 
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6.11 


6.12 


Answer: 5.23. The problem does not state whether the fault in manufacture caused the cylinder 
to distort after machining, in which case the wall thickness would be constant around the 
circumference, or both inner and outer surface were made slightly elliptical. The effect on the 
result is unlikely to be significant at such a low level of eccentricity, and for convenience the 
latter is assumed. The change in geometry is effected by modifying the original mesh to scale all 
the YNODE values by a factor of 0.99 in subprogram MESHQ immediately before finding the 
maximum dimension of the domain. Results were obtained using 4 elements per quadrant on both 
boundaries. The result for maximum hoop stress in the circular cylinder is 5 times the internal 


pressure. So, a 1% eccentricity of the cylinder gives a 4.6% increase in the maximum hoop stress. 


Answer: 299 Nm. The problem did not specify whether plane stress or plane strain is to be 
assumed. Because the state of deformation in the rubber bush is one of shearing, plane stress 
and plane strain give the same results. The problem calls for torque to be computed (via shear 
stress on the inner boundary) for a given rotation of this boundary. This rotation displacement 
boundary condition cannot be applied using the present program. On the other hand, a uniform 
shear stress can be applied to the inner boundary: a value of 1 N/m? was used. Results were 
obtained using 4 elements per quadrant on both the inner and outer circular boundaries. The 
computed value of tangential displacement at the inner boundary was 0.15750 x 107°? m. At 
a radius of 0.01 m, this corresponds to an angular displacement of 0.15750 x 107’ radians. 
The shear stress required to produce the required rotation of 0.1 radians is therefore 
shear stress = 1 X ——_—_— = 6.3492 x 10° N/m? 
This acts over a surface area of radius 0.01 m and length 0.075 m at a distance of 0.01 m from 


the axis of the shaft. Hence, the torque is 


torque = 6.3492 x 10° x 2 xm x 0.012 x 0.075 = 299 Nm 


The problem can also be solved analytically, giving a result of 299.2 Nm. 


